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The  primary  objective  of  this  research  was  to  develop 
an  accurate  mathematical  model  in  state  space  form  for  the 
Passive  and  Active  Control  of  Space  Structures  (PACOSS) 
Dynamic  Test  Article  (DTA) ,  so  that  active  control 
techniques  Ccui  be  effectively  applied  to  the  system  in  the 
future.  To  acconplish  this  goal,  accurate  system 
identification  techniques  had  to  be  found  which  would  solve 
a  multi -input  multi -output  system.  The  system 
identification  methods  chosen,  the  Eigensystem  Realization 
Algorithm  (ERA)  and  a  Backward  Least  Squares  (BLS) 
technique,  were  compared  to  a  truth  model  with  a  known  state 
space  system,  to  determine  the  accuracy  and  appliceU3ility  of 
each  method.  With  this  completed,  PACOSS  DTA  data  was 
generated  and  a  state  space  model  was  developed.  Finally, 
Bode  plots  of  the  system  model  and  the  actual  PACOSS  system 
were  compared  to  determine  the  accuracy  of  the  models 
developed.  The  BLS  method  did  not  develop  an  adequate 
model .  The  ERA  method  developed  an  accurate  plant  matrix, 
but  the  input  and  output  matrices  were  inaccurate,  resulting 
in  a  good  match  of  system  poles,  but  a  poor  match  of  system 
zeros . 
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lODELING  OP  A  LARGE  UNDAMPED  SPACE  STRUCTURE 


USING  TIME  DOMAIN  TECHNIQUES 

I.  Introduction 

Future  space  systems  are  going  to  be  much  larger  and 
require  much  more  accurate  pointing  and  vibration 
suppression  systems.  These  systems  will  push  the  state-of- 
the-art  in  control  system  technology.  But  a  control  system 
cannot  be  designed  without  an  accurate  model  of  the  system 
to  be  controlled.  This  leads  to  the  need  for  system 
identification  techniques  which  accurately  model  large  space 
structures  with  multiple,  low  frequency,  closely  spaced 
modes  (1:463) .  The  Passive  cuid  Active  Control  of  Space 
Structures  (PACOSS)  Dynamic  Test  Article  (DTA)  was  designed 
to  simulate  a  large  vindamped  space  structure  (2:1)  .  It  will 
allow  testing  of  system  identification  techniques  and  active 
control  systems,  on  an  accurate,  but  inexpensive,  earth 
based  test  bed. 

The  objective  of  this  thesis  is  to  develop  a  state 
space  model  for  the  PACOSS  DTA  to  enable  the  development  of 
active  control  systems.  This  will  be  done  using  time  domain 
techniques  and  will  result  in  a  evaluation  of  each 
technique's  effectiveness  in  developing  models  of  large 
space  structures. 
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DTA  Background 


One  of  the  most  demanding  tasks  of  a  spacecraft  control 
system  is  the  vibration  suppression  of  the  spacecraft  after 
maneuver  or  intact  with  space  debris.  In  addition,  accurate 
pointing  of  weapon  or  sensor  platforms  may  require  not  only 
micro  radian  accuracy,  but  settling  times  under  a  second. 

For  a  large  space  structure  (LSS)  with  flexible  body 
appendages,  this  will  require  extremely  advanced  control 
systems.  The  PACOSS  DTA  has  proven  to  be  an  effective  earth 
based  test  bed  to  for  researching  the  inplementation  of  both 
passive  and  active  vibration  control  techniques  for  LSS 
(2:3)  . 

The  PACOSS  DTA  was  designed  and  built  by  Martin 
Marietta  to  simulate  a  large  flexible  spacecraft  for  the 
purpose  of  control  system  design.  It  has  multiple,  closely 
spaced,  low  frequency  modes  (most  under  lO  Hertz) ,  eight 
inertial  mass  actuators  with  collocated  accelerometers,  and 
is  isolated  from  the  environment  by  eui  air  bearing  system 
(2:1) .  Physically,  it  is  shaped  like  a  primary  and 
secondary  mirror  assembly  with  two  solar  panels  (Figure  1) . 
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This  gives  it  the  shape  and  modal  characteristics  of  typical 
large  space  structures:  closely  spaced,  low  frequency 
modes,  as  well  as  low  damping  ratios. 

System  Identification  Background 

Before  control  systems  can  be  designed,  an  accurate 
model  of  the  system  must  be  developed.  Several  previous 
attempts  have  been  made  to  model  the  PACOSS  DTA.  Martin 
Marietta,  the  developers  of  PACOSS,  used  several  types  of 
finite  element  modeling  (2) .  Specifically,  they  had 
difficulty  identifying  the  modes  in  the  1  to  10  Hz  range, 
and  found  that  mode  locations  and  dairying  factors  varied  by 
as  much  as  20  percent  about  the  average  depending  upon  the 
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particular  routine  used.  These  difficulties  were  attributed 
to  the  closely  spaced  nature  of  the  model,  and  to  the 
difficulty  in  acquiring  accurate  data  on  such  a  conplex 
system  (11:1).  This  demonstrates  the  need  to  develop  new, 
and  more  accurate  system  identification  techniques  for  LSS. 

In  addition,  AFIT  students  Captain  Scott  George  used  a 
modal  analysis  technique  (7)  and  Lieutenant  Chad  Matheson 
used  a  frequency  domain  analysis  technique  (8)  to  identify 
the  system  modes.  Although  both  were  successful  at 
identifying  the  system  modes.  Captain  George  did  not  attenpt 
to  develop  a  state  space  model,  while  the  state  space  model 
developed  by  Lieutencint  Matheson  included  only  the 
collocated  input/output  pairs. 

With  both  frequency  dcxnain  and  finite  element 
techniques  already  applied  to  the  PACOSS,  one  set  of 
identification  techniques  which  have  not  been  attenpted  are 
the  time  domain  methods.  These  involve  curve  fitting  a 
transfer  function  to  a  set  of  actual  input  disturbance 
signals  sent  to  the  PACOSS  DTA  and  their  corresponding 
output  signals  received  from  the  DTA.  Specifically  the 
Bigensystem  Realization  Algorithm  (ERA)  and  the  Backward 
Least  Squares  (BLS)  method  will  be  used  to  determine  the 
system  modes,  from  which  the  plant  matrix  (or  A  matrix)  and 
input  matrix  (or  B  matrix)  will  be  formed.  A  third  method, 
called  the  Prediction  Error  Method  (PEM)  will  be  used  in 
conjunction  with  the  above  methods  to  determine  the  output 
matrix  (or  C  matrix) . 
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The  ERA  method  is  based  upon  a  singular  value 
deconposition  of  the  Hankel  matrix  formed  from  the  Markov 
parameters.  The  Markov  parameters  are  the  system  in^ulse 
responses  (6:621).  The  Backward  Least  Squares  method  uses  a 
reversed  least  squares  technique  to  separate  the  system 
modes  from  the  conputational,  or  noise,  modes  (10) . 

The  ERA  and  BLS  methods  are  specifically  designed  to 
identify  the  system  modes  and  mode  shapes  of  large  undanped 
space  structures  given  time  domain  data.  The  ERA  method  was 
used  by  Juang  and  Pappa  at  NASA  Langley  Research  Center  to 
accurately  identify  34  modes  of  the  Galileo  space  craft 
(2:624) .  Also,  Hollkanp  was  able  to  accurately  identify  10 
modes  of  the  12  meter  truss  located  at  Wright  Laboratory, 
Wright - Patterson  AFB,  Ohio,  using  both  the  ERA  and  the  BLS 
methods  (10:553) .  Although  no  attempts  were  made  to  create 
state  space  models  of  the  structures,  both  were  very 
successful  at  identifying  system  modes,  damping  factors, 
mode  shapes,  euid  modal  cutplitude  participation  factors,  on 
large  undamped  structures  similar  to  the  PACOSS  DTA.  Given 
the  system  modes  and  danping  factors,  the  plant  matrix  can 
then  be  developed. 

Finally,  PEM  is  based  upon  an  iterative  least  squares 
numerical  technique  to  model  the  system  parameters  as  well 
as  the  errors  associated  with  the  system  output  (4:75-6) .  A 
coitputer  coded  version  of  the  Prediction  Error  Method  can  be 
found  in  the  MATLAB  System  Identification  toolbox  (3:1-1) . 
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II.  Theory 


The  development  of  the  equations  used  in  this  project 
are  divided  into  four  phases.  First  is  the  develo^ent  of 
the  discrete  state  space  model  of  any  system  and  how  it 
applies  to  the  PACOSS.  This  is  followed  by  a  develo|»nent  of 
the  techniques  which  were  used  to  identify  the  specific 
state  space  model  of  the  PACOSS.  These  techniques  were  ERA, 
BLS  euid  PEM.  The  final  section  is  a  develo^xnent  of  the 
method  used  to  combine  these  techniques  in  order  to  inprove 
the  accuracy  of  the  system  model. 

Equations  of  Motion 

The  first  step  in  the  modeling  of  any  dynamic  system  is 
to  determine  the  equations  of  motion  for  that  system.  Any 
continuous  system  can  be  modeled  as  a  lumped  parameter 
spring-mass-danper  system.  For  the  PACOSS,  with  force 
inputs,  the  equation  of  motion  is: 

K5(t)  +  Gq(t)  +  lCq(t)  =  F(t)  (1) 
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where : 


X  -  Mass  Matrix  (n  x  n) 

G  -  Dancing  Matrix  (n  x  n) 

K  -  Stiffness  Matrix  (n  x  n) 
q(t)  -  displacement  vector  (n  x  i) 
q(t)  -  velocity  vector  (n  x  l) 

^(t)  -  acceleration  vector  (n  x  i) 
n  -  number  of  degrees  of  freedom 


This  can  then  be  expressed  as  the  first  order  continuous 
system: 

^(t)  =  Ax(t)  +  Bu(t) 
y(t)  =  C5c(t)  +  Du(t) 

where : 

X  -  State  Displacement  Vector  (2n  x  l) 

X  -  State  Velocity  Vector  {2n  x  l) 

Jk  -  State  TTcuisition  Matrix  (2n  x  2n) 

B  -  State  Input  Matrix  (2n  x  m) 

C  -  State  Output  Matrix  (p  x  2n) 

D  -  Feedforward  Matrix  (p  x  m) 
n  -  Number  of  Identified  Modes 
m  -  Niunber  of  Inputs 
p  -  Number  of  Outputs 


In  discrete  format,  this  equation  becomes: 

x(k  +  l)  =  Ax(k)  +  Bu(k) 
y{k)  a  Cx(k)  +  Du(k) 


(3) 
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where  k  is  the  discrete  time  increment  variable. 

The  goal  of  system  identification  when  applied  to 
control  system  design  is  to  identify  the  A,  B,  C  and  D 
matrices  which  exactly  describe  the  system  motion.  However, 
for  a  continuous  system,  the  number  of  states,  n,  is 
infinite.  So,  an  approximation  for  a  finite  value  of  n, 
must  be  made,  and  the  motion  can  therefore  only  be 
approximated.  The  required  accuracy  of  the  approximation  is 
a  key  factor  in  the  system  identification  process. 

For  a  large  flexible  space  structure  like  PACOSS,  the 
primary  concern  is  with  the  low  frequency,  non  rigid  body 
modes,  primarily  those  under  10  Hz.  So,  10  Hz  will  be  the 
cut-off  value  for  the  determination  of  the  number  of  modes. 

Blgensystem  Realization  Algorithm 

The  Eigensystem  Realization  Algorithm  was  developed  in 
1985  as  a  method  to  determine  the  modal  parameters  of  large 
flexible  multi-input  multi-output  space  structures.  The  ERA 
produces  a  minimum- order  realization  given  dynamic  test  data 
(7:620) .  This  algorithm  is  based  upon  the  singular  values 
of  the  Hankel  matrix.  To  create  the  Hankel  matrix,  the 
Markov  parameters  must  first  be  generated.  The  Markov 
parameters  are  defined  in  discrete  time  by: 

Y{k)  *  Ca'''‘b  (4) 
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The  Hankel  matrix  is: 


y(k)  y(k+ti) 

Y(ji+k)  Y(ji+k+ti) 

Y(jr-1+J^)  y(jr-l  +lt+ t^) 


Y(k+t,.i) 
y(ji  +k+t,.i) 

Y{j^.i  +k  +  t..i) 


where  ji(i=l, . . . ,r-l)  and  ti(i=l, . . . , s-l)  are  arbitrary 
integers  defined  such  that  r  is  the  number  of  modes  to  be 
determined  (both  system  modes  and  computational,  or  noised 
derived,  modes),  and  s  is  larger  then  r  (generally  3-5  times 
the  value  of  r) . 

The  first  step  to  compute  the  Hauikel  matrix  is  to 
determine  the  velocity  impulse  response  given  the  transfer 
function  generated  from  acceleration  data.  This  transfer 
function  is  defined  as: 


where : 


G.(s) 


U(s) 


(6) 


Y(s)  =  a„s“  +  +  .  .  .  +  a^s^  +  ao  (7) 

U(s)  =  bnS“  +  +  .  .  .  +  bjS^  +  bj,  (8) 


which  is  a  rational  trauisfer  function.  However,  in  state 
space  form,  a  rational  transfer  function  will  have  a  non¬ 
zero  D  matrix.  The  ERA  method  assumes  a  zero  D  matrix. 
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But,  by  integrating  the  transfer  function  (which  is 
multiplying  by  1/s  in  the  Laplace  domain) ,  the  treinsfer 
function  becomes  non -rational  and  the  D  matrix  becomes  zero. 
Numerically,  this  integration  is: 


H.(jfa>) 

jfa> 


(9) 


Then  the  in^ulse  response,  h^(t) ,  is  found  by  tcdcing  the 
inverse  Fourier  transform  of  the  frequency  response, 

Prm  this  data,  the  Markov  parameters,  and  hence  the  Hankel 
matrix,  can  be  formed. 

Once  the  full  order  Hankel  matrix  is  determined,  the 
next  step  is  to  achieve  a  minimum  realization  by  determining 
the  order  of  H„(0)  .  This  is  done  by  using  the  singular 
value  deccmposition,  defined  as: 

H„(0)  =  (10) 

where  D  is  a  diagonal  matrix  of  the  singular  values  of 
H„(0)  in  descending  order,  and  P  and  Q  are  isometric 
matrices.  The  order  of  H„(0)  is  determined  by  the  rank  of 
D.  Unfortxinately,  for  a  real  system  which  has  noise 
present,  the  singular  values  may  never  be  zero.  In  this 
case,  it  is  necessary  to  determine  some  cutoff  value  below 
which  all  the  singular  values  are  small  relative  to  the  rest 
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of  the  singular  values,  and  can  be  assiuned  zero.  This  then 
becomes  the  rank  of  H„(0)  . 

The  minimum  realization  for  the  system  is  then; 

x(k+l)  «  D^-i/3p'rH„(l)QD^-i/2x(k)+D„i/2QTEu  (H) 

y(k)  =  ETpD^i/ax(k)  (12) 

where  ET=  [Ip,  Op,  .  .  .  Op]  (6:622). 

The  final  step  is  to  determine  the  validity  of  each  of 
the  modes,  and  remove  the  noise  created  (conputational) 
modes,  leaving  only  the  actual  (system)  modes.  This  will 
result  in  a  minimum  order  system  of  purely  system  modes. 

This  order  reduction  is  done  by  defining  the  parameter 
called  the  Extended  Modal  An^litude  Coherence  (EMIVC)  .  This 
is  a  combining  of  the  controllability  ouid  observability  of 
each  mode,  presented  as  a  percent  of  completely  control leQ:>le 
and  observable.  The  EMAC  is  a  function  of  the  singular 
values,  the  modal  participation  factors  and  the  mode  shapes 
(6:623)  and  is  used  to  help  separate  the  computational  modes 
frcmi  the  system  modes.  The  primary  assumption  behind  the 
EM^C  is  that  only  system  modes  will  be  controllable  and 
observed>le.  Therefore,  the  computational  modes  will  have 
significantly  lower  values  for  the  EMAC.  This  allows  for 
quick  and  easy  separation  of  the  system  euid  computational 
modes.  Additionally,  each  system  mode  has  a  different  level 
of  controllability  and  observability,  resulting  in  different 
value  of  the  EMAC  for  each  mode.  Given  a  limit  to  the 
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nxunber  of  syst^  modes  which  can  be  kept,  this  provides  a 
qucuititative  measurement  to  determine  which  modes  should  be 
kept  and  which  modes  can  be  discarded. 

Backwards  Least  Squares 

The  second  method  used  for  system  identification  was  a 
multi- input,  multi -output  Backwards  Least  Squares  technique 
developed  by  Hollkair^)  (10:549).  The  forward  least  squares 
method  is  based  upon  the  fact  that  a  time- series  model  can 
be  expressed  in  the  form: 

y(k)  =  -Ajy(k-l) -Aay{k-2)  -  . .  . -A^(k-p) 

+B,u(k)  +B,^u(k-1)  +  . .  .+BpU(k-p)  (13) 

where  y(k)  is  the  k^‘‘  output  response,  u(k)  is  the  ktJ>  input, 
p  is  the  order  of  the  model,  the  A^,  A^,  ...,  A^  and  B^,  B,, 

.  .  . ,  Bq  matrices  are  the  matrices  of  model  parameters 
(10:549).  Then,  for  an  over  determined  system,  the  solution 
is: 


Y  =  (14) 

where  Y  is  the  vector  of  output  data,  F  is  a  matrix  of 
output  and  input  data,  and  B  is  the  vector  of  unknown 
parameters  (10:550) .  For  the  forward  least  squares 
technique,  the  solution  to  this  problem  will  be  a  system 
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model,  but  it  will  not  be  minimum  order.  That  is  because  6 
will  include  both  system  (actual)  and  con^utational  (error 
created)  modes.  Order  reduction  of  the  system  can  be 
sin^lified  by  performance  of  the  Backward  Least  Squares 
technique  (10:550) . 

Any  time- series  system  can  be  written  in  a  backwards 
manner  of  the  form: 

y(k)  =  -Ciy(k+1) -Cay(k+2)  - .  .  . -C^(k+p) 

+DjU(k) +DjU(k+l)  + . .  .+DpU(k+p)  (15) 

where  the  C^,  Cj,  .  .  . ,  and  D^,  Dj,  ...»  matrices  are  the 
bac)cward  parameter  matrices.  The  solution  for  the  unknown 
parameters,  0,  is  found  the  same  as  for  the  forward  method, 
using  equation  14.  Unlike  the  forward  least  squares  method, 
the  backward  method  will  force  all  the  computational 
eigenvalues  inside  the  unit -amplitude  circle  (10:550) .  This 
means  that  only  the  eigenvalues  outside  the  unit -amplitude 
circle  need  be  kept,  as  they  are  the  system  modes.  In  this 
way,  a  reduced  order  system  is  found.  The  conversion  back 
to  forward  least  squares  (to  create  a  stable  system  with 
eigenvalues  less  than  one)  is: 

K  =  (16) 

Bi  =  Cp-"Dp-i  (17) 
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The  new,  reduced  order  suid  matrices  Ccui  now  be  used  to 
solve  for  the  new  reduce  order  system  modes. 


Prediction  Error  Method 

The  ERA  and  BLS  methods  were  developed  primarily  to 
develop  natural  frequencies,  dan5>ing  factors,  mode  shape 
matrices  and  modal  participation  matrices,  (e.g.  a  modal 
model  of  the  system)  not  a  state  space  model.  Therefore,  an 
additional  technique  had  to  be  added  to  enable  development 
of  cui  accurate  system  in  state  space  form.  If  one  assumes 
that  the  ERA  and  BLS  methods  are  able  to  obtain  accurate  A 
and  B  matrices,  then  the  PEM  has  the  ability  to  ccwnpute  the 
C  matrix  given  the  A  cuid  B  matrices. 

The  Prediction  Error  Method  is  found  in  the  system 
identification  toolbox  created  by  Math  works  Inc.  for  MATLAB 
(3)  .  This  method  atten?)ts  to  minimize  the  error  gradient 
between  the  model  generated  output  and  the  actual  system 
output.  Given  a  set  of  input -output  data,  and  a  model 
format,  with  a  set  of  fixed  parameters  and  a  set  of  unknown 
parameters,  this  program  attempts  to  vary  the  unknown 
parameters  along  the  gradient  direction  to  fit  the  model 
results  to  the  known  set  of  data.  The  accuracy  and 
confutation  time  for  this  routine  vary  with  model  size  and 
the  number  of  unknown  parameters. 
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This  technique  uses  an  iterative  Gauss -Newton  algorithm 
to  solve  the  equation: 

^Mi(q)y(t)  =  |-^u(t  -  nk)  e(t)  (18) 

F„(q)  D«(q) 

Where  e(t)  is  the  measurement  error,  u{t-nk)  is  the  input 
data,  y(t)  is  the  output  data  and  ApB^^q)  <  ^paf(q) «  Cpai(q)  > 
DpB,(q)  ,  and  PpBM(q)  are  of  the  form: 

Xpa,(q)=I+aiq-i+.  .  .+a,q-“  (19) 

This  method  solves  for  a  particular 

d=[ai,  . .  .a„,bi,  .  .  .b„,Ci, . .  .Ca,di,  . .  .d^/fi,  •  •  -fj  (20) 

which  minimizes  the  square  of  the  error  between  the  output, 
y,(t),  of  the  actual  system  and  the  output,  y„(t)  ,  of  the 
model  (3:1-7)  using  a  gradient  minimization  technique.  The 
actual  system  is  assumed  be  a  constant  coefficient  system. 

Method  Comparison  and  Combination 

The  BRA  and  the  BLS  methods  were  the  primary  techniques 
used  to  confute  the  system  parameters.  In  addition,  since 
they  are  different  methods  which  use  different  types  of  data 
to  generate  the  same  solution,  their  results  were  compared 
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to  help  determine  the  validity  of  the  final  system  model 
developed. 

The  plant  matrix  was  created  by  generating  the  discrete 
eigenvalues  using  the  equation: 

z  =  exp((-i'w„  +  jWd)/f,]  (21) 

where  z  is  the  eigenvalue,  f  is  the  damping  factor,  is 
the  undanped  natural  frequency,  is  the  danped  natural 
frequency,  and  f,  is  the  sampling  rate  (2:622).  The 
discrete  eigenvalues  then  become  the  diagonals  of  the  plant 
matrix  in  discrete  time.  Since  this  is  a  conplex  matrix,  it 
must  be  transformed  into  the  real  2x2  block  diagonal  form. 
This  is  then  the  real  A  matrix  which  can  be  used  in  control 
system  design.  The  cotrplex  modal  participation  matrix  is 
the  continuous  time  form  of  the  B  matrix.  This  B  matrix 
must  be  transformed  to  real,  discrete  time  form,  to  be 
consistent  with  the  A  matrix. 

That  leaves  only  the  C  matrix  as  an  unknoim.  PEM  was 
used  to  compute  an  accurate  C  matrix.  This  technicjue 
beccmies  more  effective  the  fewer  the  parameters  there  are  to 
estimate,  since  PEM  computes  a  gradient  for  each  unknown 
parameter  emd  searches  along  that  direction  in  an  attempt  at 
minimization.  So  with  the  A  and  B  matrices  as  known 
parameters,  all  that  PEM  needs  to  approximate  is  the  C 
matrix.  By  approximating  the  PACX)SS  as  a  linear  system, 
then  each  output  is  independent  of  the  other  outputs.  Then, 
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PACOSS  can  be  approximated  as  multiple  single  output 
systems.  The  number  of  parameters  PEM  needs  to  estimate  is 
then  equal  to  a  single  row  of  the  C  matrix.  So,  PEM  was 
performed  8  times  for  each  of  the  8  outputs,  then  recombined 
into  a  single  C  matrix  of  the  form: 


This  reduced  the  niomber  of  parameters  PEM  needed  to 
estimate,  while  providing  an  accurate  C  matrix  of  the 
system.  The  combined  methods  of  ERA  A  and  B  determination 
and  PEM  C  determination  will  be  termed  ERA/PEM  in  the 
remainder  of  the  thesis.  The  combined  method  of  BLS  A  and  B 
determination  and  PEM  C  determination  will  be  termed  BLS/PEM 
in  the  remainder  of  the  thesis. 
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Before  attenpting  to  identify  the  PACOSS  DTA,  the 
accuracy  of  the  identification  methods  selected  were 
validated  on  a  known  truth  model.  This  required  a  spring- 
mass  -danger  system  with  ]aiown  plant,  input  euid  output  to  be 
created.  ERA/PEM  cuid  BLS/PEM  were  then  applied  to  this 
known  system  amd  the  models  developed  were  conpared  to  the 
truth  system  to  determine  the  accuracy  of  each  method. 


To  determine  the  accuracy  and  applicability  of  the 
system  identification  techniques  selected,  it  was  decided  to 
apply  them  to  a  low  order  "truth"  model  first.  This  would 
allow  a  conparison  of  each  model  against  known  A,  B,  C,  and 
D  matrices.  To  accurately  determine  the  applicability  of 
the  system  identification  techniques  selected,  the  truth 
model  must  contain  the  same  characteristics  as  the  system  to 
be  identified,  in  this  case,  the  PACOSS  DTA.  These 
characteristics  include  multi -input  multi-output  and 
closely  spaced  poles  at  or  below  10  Hz  with  low  dancing 
(2:15).  As  a  result,  a  spring-mass -dan5)er  system  of  the 
form  of  equation  (1)  was  developed  (Figure  2) .  This  was  a  2 
input,  2  output  system  with  masses  of: 
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ml  =  0.61743  kg 
in2  =  0.19957  kg 
m3  s  1.6582  kg 

dancing  of: 

cl  =  0.031633  kg/sec 
C2  =  0.32738  kg/sec 
c3  =  0.32918  kg/sec 
c4  =  0.019153  kg/sec 


and  spring  consteuits  of: 

kl  =  52.058  N/m 
k2  =  0.89183  N/m 
k3  =  19.756  N/m 
k4  =  178.08  N/m 


Figure  2:  Truth  Model 
Physical  System 


The  outputs  were  position  measurements  of  each  mass. 
They  were  combined  into  discrete  time  state  space  format 


19 


(Equation  3)  .  The  natural  frequency  euid  danping  factors  for 
a  sanpling  rate  of  14.2  samples  per  second  are  presented  in 
Table  1. 


Table  1:  Truth  Model  Modes 


(Oa  (Hz) 

r  (%) 

1.4414 

1.4537 

1.8817 

2.0162 

9.9314 

8.0431 

Prom  which  the  following  discrete  time,  state  space  system 
was  developed: 


0.7999 

0.0005 

0.0026 

0.0639 

0.0012 

0.0001 

0.0022 

0.7766 

0.1950 

0.0036 

0.0575 

0.0084 

0 

0.0255 

0.7250 

0 

0.0010 

0.0630 

-5.4780 

-0.0279 

0.1048 

0.7647 

0.0305 

0.0045 

-0.0542 

-5.8419 

4.6834 

0.0944 

0.5910 

0.2881 

0.0015 

0.6453 

-7.4133 

0.0017 

0.0347 

0.7134 

0.0001 

0 

0.0109 

0.0001 

0.0001 

0.0014 

0.0059 

0.0001 

0.2881 

0.0051 

0.0051 

0.0.380 

0  10 

0  0  0“^ 

0  0  1 

0  0  0_ 
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Which  are  the  A,  B,  and  C  matrices  defined  in  Eq.  (3) .  The 
D  matrix  is  zero. 


Truth  Model  Identification 

The  identification  of  a  known  truth  model  was  the  first 
step  in  euialyzing  the  effectiveness  of  the  system 
identification  techniques  chosen.  For  the  ERA  method,  an 
impulse  response  of  the  truth  model  was  generated.  A  zero 
mean  Gaussleui  measurement  error  (RMS  value  of  O.i)  was  added 
to  the  output.  The  ERA  method  used  a  20  column  (lO  modes) 
and  60  row  (58  data  points  for  4.06  seconds  of  data)  Hankel 
matrix.  Only  the  18  largest  singular  values  of  the  Hankel 
matrix  were  retained,  with  the  two  remaining  singular  values 
discarded  due  to  there  relatively  small  value. 

In  a  similar  manner,  the  BLS  method  had  a  random  signal 
of  zero  mean  Gaussieui  noise  which  had  an  RMS  value  of  l, 
used  as  an  input  to  the  truth  model.  A  measurement  error  of 
zero  mean  Gaussian  noise  with  an  RMS  value  of  O.Ol,  was 
added  to  the  output.  The  resulting  natural  frequencies  and 
danqping  factors  confuted  by  each  method  are  shown  in  Table  2 
along  with  the  actual  natural  frequencies  emd  damping 
factors  of  the  system. 

From  these  results,  it  is  evident  that  the  ERA  method 
was  able  to  Identify  the  actual  modes  of  the  system  with  a 
great  deal  of  accuracy,  despite  having  fairly  high 
measurement  errors  (approximately  10%  of  the  input) .  It  was 
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even  able  to  separate  and  Identify  the  closely  spaced  modes. 


Table  2:  Identified  Truth  Model  Modes 


Actual 

ERA 

BLS 

(On  (Hz) 

r  {%) 

(On  (Hz) 

r  (%) 

(On  (Hz) 

f  (%) 

1.4414 

2.0162 

1.4414 

2.0162 

1.4063 

6.4529 

1.4537 

9.9314 

1.4537 

9.9314 

- 

- 

1.8817 

8.0431 

1.8817 

8.0431 

1.9066 

3.4481 

The  BLS  technique  did  not  identify  the  modes  or  dan^jing 
factors  as  accurately  as  the  ERA  technique,  despite  having 
considerably  lower  measurement  error  noise.  This  is 
discouraging  since  the  PACOSS  DTA  has  several  closely  spaced 
poles,  and  a  considerable  amount  of  noise. 

Next,  the  A  and  B  matrices  generated  by  each  method 
were  inserted  into  PEM  to  compute  the  C  matrix.  This  gave  a 
complete  state  space  model  of  the  truth  system  (Appendix  A) . 
Now,  the  time  and  frequency  response  outputs  for  each  model 
were  compared  to  those  generated  by  the  truth  model 
(Appendix  A) .  Figure  3  shows  that  the  ERA/PEM  method 
produces  a  virtually  identical  frequency  response  (dashed 
line)  to  the  truth  system  (solid  line)  below  5  Hz,  with  some 
error  entering  into  the  system  between  5  and  10  Hz.  What  is 
inqportant  to  note  is  the  accuracy  with  which  the  poles  and 
zero  were  identified. 
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Frequency  (Hz) 


Figure  3:  ERA/PEM  Magnitude  Plot 

It  is  also  iir^ortant  to  note  that  the  ERA  method  without  PEM 
did  not  identify  the  zeros  of  the  system  correctly.  As  can 
be  seen  in  Figure  4,  the  ERA  method  alone  added  a  zero  not 
found  in  the  actual  system.  The  PEM  routine  was  needed  to 
accurately  identify  the  system  zeros. 


Frequency  (Hz) 


Figure  4:  ERA  Magnitude  Plot 


The  results  of  the  BLS/PEM  method  were  not  as 
encouraging.  By  examining  the  frequency  response  plots,  it 
is  evident  that  the  BLS/PEM  model  identified  two  closely 
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spaced  poles  as  a  single  pole  (Figure  5) .  The  poles  and 
zeros  of  the  BLS  model  (dashed  line)  did  not  match  up  well 
with  those  of  the  true  system  (solid  line) ,  and  the  BLS 
model  added  a  zero  at  6  Hz,  and  failed  to  identify  the  pole 
at  1.4537  Hz.  The  resulting  magnitude  plot  shows 
significant  error.  This  was  with  only  1%  measurement  error. 


Appendix  A  contains  the  state  space  model  developed  by 
both  ERA/PEM  and  BLS/PEM,  as  well  as  all  8  bode  plots. 

Once  the  system  identification  techniques  had  been 
tested  on  the  truth  model,  the  next  step  was  to  use  these 
techniques  to  confute  a  state  space  model  for  the  PACOSS 
DTA. 
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IV.  Experimental  Procedure 


The  identification  of  the  system  model  for  the  PACOSS 
DTA  was  divided  into  three  phases.  The  first  phase  was  to 
evaluate  different  time  domain  system  identification 
techniques  to  determine  which  was  the  most  appropriate.  The 
results  of  phase  one  are  detailed  in  Chapter  III.  The 
second  phase  was  to  actually  acquire  the  input  and  output 
data  from  the  PACOSS  system.  While  the  third  phase  was  to 
run  the  actual  data  through  the  chosen  system  identification 
techniques  to  determine  the  system  model  for  the  PACOSS  DTA. 

Data  Generation 

Once  the  accuracy  of  each  identification  technique  was 
determined,  the  actual  PACOSS  input/output  data  was  then 
generated.  First,  all  equipment  (in  particular  the 
accelercaneters)  had  to  be  calibrated  to  identify  emy  bias 
allowing  the  generation  of  accurate  and  consistent  data.  A 
Techtronix  Digital  Analyzer  2642  was  used  to  generate  the 
random  input  signals  amd  to  record  the  output  from  the 
PACOSS.  It  also  computed  and  recorded  frequency  response 
functions  for  individual  input/output  pairs  of  the  PACOSS 
DTA.  One  limitation  of  the  Digital  Analyzer  is  that  it  can 
only  simultaneously  record  2  channels  of  data.  The  PACOSS 
system  is  8  input  and  8  output,  therefore  this  causes  some 
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difficulty  in  acquiring  the  data.  Since  the  BLS  cuid  PEM 
techniques  require  knowing  the  input  corresponding  to  each 
output,  1  input  and  l  output  (velocity)  channel  were 
recorded,  requiring  64  separate  data  runs.  This  allowed  the 
combination  of  all  64  data  samples  to  be  combined  into  a 
single  8  input  8  output  model. 

The  next  critical  decision  was  to  determine  the 
sair5)ling  rate.  San^ling  rate  is  in?>ortant  because  it 
determines  the  maximum  frequency  of  the  system  which  can  be 
identified  (4:378).  In  addition,  a  high  sanpling  rate  is 
necessary  to  distinguish  between  closely  spaced  poles. 
However,  it  is  in5)ortauit  to  remember  that  the  higher  the 
sampling  rate,  the  more  data  there  is  generated.  This  is  an 
important  trade-off  since  there  are  64  input/output 
combinations  for  which  data  must  be  generated.  In  order  to 
generate  sufficient  data  yet  maintain  a  Nyquist  frequency  in 
the  10  Hz  region,  the  random  input  response  data  was 
acquired  at  25.6  Hz  for  20  seconds  of  data.  The  frequency 
response  function  plots  were  generated  for  data  from  0  to  10 
Hz,  using  the  Sweptsine  sub  function  on  the  2642  Analyzer 
with  a  step  size  of  0.025  Hz  collecting  401  points  of  data. 
This  data  was  collected  by  Lieutenant  Chad  Matheson  as  part 
of  his  thesis,  in  1992  (8) . 

Since  the  data  collected  by  Lieutencuit  Matheson  was  a 
set  of  frequency  response  functions  created  frcati 
acceleration  data,  and  the  ERA  method  requires  impulse 
response  data,  with  output  measurements  of  velocity,  not 
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acceleration,  a  trcuisformatlon  of  the  data  had  to  be  made. 
The  first  step  was  to  integrate  the  frequency  response 
function  data.  In  the  frequency  domain,  this  becomes: 


Hv(ja)) 


H,(ja>) 

joj 


(23) 


The  primary  problem  encountered  was  to  determine  Hv(ja))  at 
w  =  0  Hz  (DC) .  Since  most  of  the  data  below  l  Hz  was  not 
very  accurate,  and  the  modes  below  l  Hz  include  primarily 
rigid  body  and  pendulum  modes  which  will  not  be  controlled, 
it  was  decided  to  make  Hv(0)  =  0.  This  is  the  same  method 
used  by  the  MATH  subroutine  on  the  Techtronix  Digital 
Analyzer  2642. 

Once  the  frequency  response  functions  with  respect  to 
velocity  were  created,  the  inverse  fourier  transform  of  the 
data,  could  then  be  taken,  again  using  a  MATLAB  script  file 
(see  Appendix  D) .  This  script  file  used  the  IPFT  function 
in  MATLAB  for  2048  points,  emd  multiplied  the  results  by  the 
time  scaling  factor,  which  is  2048.  This  then  provided  the 
velocity  impulse  response  which  the  ERA  method  requires. 

The  other  set  of  data  required  was  the  random  input 
data.  This  data  was  generated  by  sending  a  random  input 
signal  to  the  proof  mass  actuators  to  excite  the  DTA.  The 
acceleration  output  signal  thus  generated  was  sent  through 
the  analog  integrator  located  in  the  motor  control  units, 
and  was  recorded  by  the  Analyzer.  However,  only  the 
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collocated  velocity  data  was  valid  since  the  velocity  signal 
for  the  non- collocated  measurements  was  so  small  as  to  be 
indistinguishable  from  the  ambient  noise  of  the  system. 

This  was  true  only  for  velocity  data  gathered  from  the 
Analyzer.  This  is  because  although  the  accelerations  are 
large,  they  occur  for  only  a  short  time  in  each  direction, 
resulting  in  small  velocities.  The  accelerations  for  the 
collocated  data  are  an  order  of  magnitude  larger  than  for 
non- collocated  data,  resulting  in  collocated  velocity  data 
which  is  distinguishcdsle  from  the  ambient  noise  of  the 
system.  For  this  reason,  the  data  inserted  into  PEM  to 
generate  the  C  matrix  only  included  the  collocated  rate 
data. 

Model  Generation 

with  the  data  generated,  the  next  step  is  to  identify 
the  system  model .  The  in5)ulse  response  data  was  inserted 
into  the  ERA  software  while  one  set  of  random  input 
responses  (and  corresponding  input  signals)  were  inserted 
into  the  BLS  software.  Both  methods  generated  natural 
frequencies,  dancing  factors,  mode  shape  matrices  and  modal 
participation  matrices.  The  plant  matrix  was  generated  from 
the  natural  frequencies  and  dan^jing  factors  while  the  input 
matrix  was  generated  by  taking  the  discrete  transformation 
of  the  modal  participation  matrix.  These  matrices  were  then 
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transformed  to  real  matrices  (2x2  block  diagonal  form  for 
the  A  matrix)  . 

The  only  matrix  left  to  identify  was  the  output  matrix. 
This  was  were  the  PEM  method  was  employed,  /mother  set  of 
random  input  response  data  was  inserted  into  the  PEM  routine 
(see  Appendix  D) .  The  elements  of  the  A  and  B  matrices  were 
fixed,  while  the  PEM  routine  was  allowed  to  determine  the 
elements  of  each  row  of  the  C  matrix  (given  the  real  form  of 
the  mode  C  matrix  derived  by  the  ERA  or  BLS  procedure  as  an 
initial  condition) .  By  fixing  the  A  and  B  matrices,  the  PEM 
routine  was  able  to  accurately  determine  the  C  matrix 
because  of  the  relatively  few  numbers  of  variable  parameters 
it  had  to  determine. 

The  final  step  in  the  system  identification  process  was 
to  test  the  accuracy  of  the  model.  Since  the  numerous 
errors  associated  with  time  domain  data  make  conparison  in 
the  time  domain  difficult,  and  because  small  errors  in  the 
frequency  domain  (such  as  not  identifying  one  pole  or  zero) 
translate  to  large  errors  in  the  time  domain,  it  was  decided 
that  only  a  comparison  of  the  system  transfer  functions,  and 
no  conqparison  of  time  domain  output  would  be  made.  A 
comparison  of  the  frequency  response  of  each  model  was 
con^ared  to  the  frequency  response  functions  of  the  DTA. 

The  PACOSS  system  frequency  response  functions  were  the  same 
velocity  transfer  functions  which  were  used  earlier  in  the 
generation  of  the  impulse  response  data  for  the  ERA.  The 
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results  of  the  PACOSS  identification  are  presented  in 
Chapter  V. 

With  this  identification  accomplished,  there  is  a 
plant,  Input,  and  an  output  matrix  describing  the  motion  of 
the  PACOSS  DTA.  This  state  space  model  can  be  used  to 
design  an  active  control  system  for  the  PACOSS. 
Additionally,  there  is  a  conparison  of  two  system 
identification  techniques  for  their  accuracy  in  modeling 
large,  multi- input  multi-output,  undanped  space  structures 
with  multiple,  closely  spaced,  low  frequency  modes. 
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with  the  model  identification  techniques  selected,  and 
a  validation  of  their  ability  to  generate  accurate  state 
space  models,  the  final  step  is  to  generate  a  state  space 
model  for  the  PACOSS  DTA.  Both  the  ERA/PEM  and  BLS/PEM 
methods  were  used,  and  the  results  were  compared  to  the 
frequency  response  of  the  PACOSS  DTA  to  determine  the 
accuracy  of  model  thus  developed. 

BRA  Identification 

Since  the  ERA  method  produced  an  extremely  accurate 
model  of  the  truth  system,  it  was  the  first  method  employed 
to  identify  the  PACOSS  DTA.  The  impulse  response  data  was 
inserted  into  the  ERA/PEM  software,  using  a  Hankel  matrix 
that  had  90  columns  (45  modes)  and  400  rows  (79  data  points 
for  1.66  second  of  data) .  The  singular  value  cutoff  was  at 
the  86*^*'  singular  value.  The  A  and  B  matrices  thus  created 
were  inserted  into  PEM  and  a  C  matrix  was  conputed.  Appendix 
B  contains  the  state  space  model,  in  discrete  time  form, 
developed  by  ERA/PEM. 

As  can  be  seen  in  TeUDle  3  ERA/PEM  identified  lO  modes 
under  the  lO  Hz  cutoff.  These  modes  were  closely  space  and 
they  all  had  low  damping  ratios. 
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Table  3:  ERA/PEM  Identified 
PACOSS  Modes 


ERA/PEM 

(On  (HZ) 

r  (%) 

1.5258 

5.1503 

2.8969 

2.0358 

3.9228 

1.2975 

4.7210 

1.7325 

4.9421 

0.9926 

6.9764 

5.3497 

7.2674 

0.4309 

9.1699 

0.4153 

9.2688 

4.0237 

9.6443 

1.8441 

In  order  to  test  the  accuracy  of  these  modes,  a 
conparison  was  made  to  the  experimentally  obtained  frequency 
response  function  plots  from  the  PACOSS.  Figure  6  is  a 
comparison  of  the  PACOSS  frequency  response  plots  (solid 
line)  and  the  model  developed  using  the  C  matrix  generated 
by  ERA  only  (dashed  line) .  It  can  be  noted  that  few  zeros 
were  identified  and  many  of  the  poles  were  ceuiceled  out  by 
the  zeros  that  were  identified. 
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Frequency  (Hz) 

Figure  6:  ERA  Magnitude  Plot  for 
Actuator  l/Accelerraneter  l 


Figure  7,  Is  a  ccxiparlson  of  the  PACX)SS  frequency  response 
plots  (solid  line)  and  the  model  which  results  from  the  PEM 
developed  C  matrix  (dashed  line) .  This  Is  a  more  accurate 
match  with  no  poles  canceled  by  zeros.  Appendix  C  contains 
the  graphs  of  all  the  actuator/accelerometer  pair 
conparlsons  between  each  PACX>SS  frequency  response  plot  and 
Its  corresponding  ERA/PEM  model  frequency  response  plot. 


Frequency  (Hz) 

Figure  7:  ERA/PEM  Magnitude  Plot  for 
Actuator  l/Accelerometer  l 
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The  ERA/PEM  method  shows  a  good  correlation 
between  l  Hz  and  10  Hz  for  the  system  poles,  but  not  for  the 
system  zeros.  The  ERA/PEM  method  identified  almost  all  the 
poles,  and  many  of  the  zeros.  Which  specific  poles  cuid 
zeros  it  identifies  accurately  varies  with  the 
actuator /accelerometer  pair  being  examined.  As  can  be  seen 
in  Figure  7,  for  the  actuator  1/accelerometer  l  pair,  most 
of  the  poles  are  identified  by  ERA/PEM  quite  accurately,  but 
4  significant  zeros  are  not  identified.  This  is  typical  of 
the  system  model  developed  (see  Appendix  C  for  other 
actuator /accelerometer  pairs) . 

What  this  signifies  is  that  the  A  matrix  for  the  system 
model  is  accurate  and  could  be  used  in  control  system 
design.  The  B  and  C  matrices,  however,  are  not  accurate. 
They  do  identify  some  zeros,  but  they  miss  several 
significant  zeros.  Overall,  the  ERA/PEM  method  does  a  good 
job  of  modeling  a  very  complex  system,  but  the  model 
developed  would  not  be  accurate  enough  to  enable  control 
system  design. 

BLS  Identification; 

Although  the  BLS  method  was  not  as  accurate  as  the  ERA 
method  in  identifying  the  truth  model,  in  particular  when 
there  was  a  large  measurement  noise,  it  was  thought  that  BLS 
could  help  to  confirm  some  of  system  modes  developed  by  ERA. 
For  this  method,  the  set  of  random  input  data  was  used  to 
develop  the  system  model.  As  before,  the  A  and  B  matrices 
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were  inserted  Into  PEM  to  confute  a  C  matrix  and  provide  a 
complete  state  space  model  of  the  PACOSS  DTA. 

As  can  be  seen  in  Table  4,  the  BLS  method  identified  10 
modes.  These  modes  were  grouped  aro\ind  1.5  Hz  and  12  Hz, 
with  most  having  danping  ratios  in  excess  of  30%. 


Table  4:  BLS/PEM  Identified 
PACOSS  DTA  Modes 


BLS/PEM 

Cdn  (HZ) 

f  (%) 

0.7202 

81.3024 

1.2406 

58.2258 

1.5165 

62.8178 

1.4941 

59.9581 

1.5001 

53.0558 

1.3518 

30.8457 

1.4993 

40.9521 

12.6124 

3.5030 

12.6761 

1.6485 

12.8435 

8.2277 

These  modes  are  very  different  from  the  actual  system 
response  as  seen  in  Figure  8,  where  the  modes  are  spread  out 
from  2  Hz  to  10  Hz  with  danping  ratios  under  5%.  As  can  be 
seen,  there  is  almost  no  correlation  between  the  BLS/PEM 
frequency  response  (dashed  line)  and  the  actual  PACOSS  DTA 
frequency  response  (solid  line) . 
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Frequency  (Hz) 


Figure  8:  BLS  Magnitude  Plot  for 
Actuator  #i/Accelerometer  #i 

This  lack  of  correlation  leads  to  the  conclusion  that 
the  state  space  model  developed  using  the  BLS/PEM  method  was 
not  accurate  enough  to  enable  the  development  of  a  control 
system. 

Sources  of  Error 

The  primary  reason  the  model  did  not  exactly  match  the 
PACOSS  DTA  in  the  frequency  domain  was  because  of  the 
numerous  errors  inherent  in  the  system  identification 
process,  in  addition  to  the  errors  particular  to  this 
specific  experiment.  Inherent  in  the  attenpt  to  identify 
any  real  world  system  is  the  fact  that  any  real  world  system 
is  a  continuous  system  with  eui  infinite  number  of  modes. 
However,  a  model  is  approximated  by  a  system  of  discrete 
point  masses  euid  can  only  have  a  finite  nxunber  of  modes. 

The  atten^t  is  made  to  retain  the  modes  which  have  the  most 
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effect  on  the  system  behavior.  But,  the  fewer  modes 
retained,  the  less  accurate  the  model.  The  ERA/PEM  method 
generated  a  model  with  10  modes.  This  is  far  from  an 
infinite  number  of  modes,  but  most  of  the  power  is  located 
in  these  modes. 

Another  set  of  errors  inherent  in  any  system 
identification  process  is  the  noise  introduced  into  the 
response  measurements.  This  noise  can  come  from  a  number  of 
sources  and  includes  accelerometer  noise  and  bias,  errors 
measuring  the  inputs,  air  currents,  modes  introduced  by  the 
attachment  points  between  the  experiment  and  the 
environment,  and  rigid  body  or  pendulum  modes.  These  noises 
all  very  in  magnitude,  and  are  not  always  zero  mean 
Gaussian.  This  makes  them  difficult  to  predict  and  remove 
from  the  data  before  processing.  Since  most  system 
identification  models  attenpt  to  account  for  zero  mean 
Gaussian  noise,  noise  which  does  not  fit  this  description 
will  effect  the  system  model.  The  result  is  that  system 
identification  techniques  model  this  non- zero  mean  non- 
Gaussian  noise  as  additional  system  modes,  therefore 
creating  a  model  with  more  modes  than  the  actual  system. 

For  this  particular  experiment,  there  were  several 
known  sources  of  noise.  The  first,  and  probably  the  most 
significant,  was  accelerometer  noise  and  bias.  For  velocity 
data,  the  accelerometer  ambient  noise  level  was  on  the  order 
of  10  mVolts,  with  a  steady  state  bias  ranging  from  -lOO 
mVolts  to  -400  mVolts.  This  was  significant  because  the  non 
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collocated  velocity  data  had  velocity  values  on  the  order  of 
10  mVolts.  This  meemt  that  none  of  the  non  collocated 
velocity  measurements  could  be  used  because  the  data  could 
not  be  distinguished  from  the  noise.  Additionally,  the 
PACOSS  air  bearing  isolation  system  has  a  known  natural 
frequency  of  0.6  Hz,  which  was  accounted  for,  and  an  unknown 
nvimber  of  additional  frequencies  which  introduce  noise  into 
the  output  response.  Additionally,  there  were  low  frequency 
pendulum  modes  which  effected  the  data.  For  this  reason, 
the  data  below  i  Hz  was  viewed  with  skepticism. 
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YJL. _ Conclusions  and  RecommendatioriB 


The  final  result  of  this  research  is  that  the  system 
identification  methods  used  were  unable  to  create  a  state 
space  model  with  sufficient  accuracy  to  enable  the  design  of 
an  active  control  system  for  the  PACOSS  DTA.  They  were, 
however,  cUole  to  accurately  identify  the  system  modes  and 
damping  factors,  from  which  the  plant  matrix  was  created. 

Conclusions 

The  ERA  method  appears  to  be  a  accurate  method  for 
determining  system  modes.  It  accurately  predicted  the 
system  poles,  as  is  evident  in  the  frequency  domain  plots  in 
Appendix  C.  This  is,  after  all,  what  it  was  designed  to  do. 
One  reason  this  method  had  difficulty  identifying  all  the 
modes  at  all  the  sensor  locations  was  due  the  difficulty  it 
had  in  identifying  the  system  zeros.  For  the  most  part, 
BRA/PEM  failed  to  identify  the  system  zeros,  resulting  in  a 
state  space  model  with  far  fewer  zeros  than  were  foiind  in 
the  actual  system.  The  result  of  this  may  be  that  at 
certain  actuator /accelerometer  pairs,  the  ERA/PEM  method  did 
not  identify  a  zero,  which  in  the  actual  system  canceled  a 
pole.  The  model  frequency  response  would  then  show  a  pole 
where  the  PACOSS  frequency  response  does  not,  due  to  pole- 
zero  cancellation. 
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The  inability  of  the  ERA/PEM  method  to  accurately 
identify  all  the  system  zeros  was  disappointing,  but  not 
unexpected,  since  the  method  was  designed  as  a  modal 
identification  technique,  not  a  state  space  model 
development  technique.  But,  for  control  system  design,  the 
location  of  the  zeros  is  as  i:..portant  as  the  location  of  the 
poles.  For  this  reason,  to  develop  a  truly  accurate  active 
control  system  for  the  PACOSS,  some  other  method  will  have 
to  be  used,  instead  of,  or  in  conjunction  with,  the  ERA 
method . 

This  thesis  atten^jted  to  supplement  the  ERA'S  zeros 
finding  capabilities  with  the  PEM  method.  For  the  small 
truth  model,  this  worked  well.  But  for  the  larger,  amd  more 
ccm^lex  PACX)SS  system,  this  method  did  not  work  very  well. 
PEM  did  identify  several  zeros,  but  it  missed  several 
others.  Poor  performance  in  the  larger  system  is  most 
likely  due  to  the  method  PEM  employs  in  finding  zeros.  It 
uses  a  gradient  search  to  identify  unknown  parameters  of  the 
A,  B,  and/or  C  matrix.  For  this  thesis,  the  parameters  of 
the  C  matrix  were  identified  one  row  at  a  time.  This  meams 
PEM  was  trying  to  minimize  the  gradient  for  20  unJmown 
parameters,  as  opposed  to  the  6  unknown  parameters  in  the 
truth  model. 

To  help  verify  that  this  set  of  zeros  was  the  best 
result  PEM  could  achieve,  the  C  matrix  was  held  constant  amd 
the  B  matrix  was  identified  one  column  at  a  time.  The 
results  were  that  the  same  set  of  zeros  were  identified  amd 
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were  not  identified  as  before.  Finally,  when  both  B  and  C 
were  inserted  into  PEM  and  allowed  to  vary,  the  results  were 
worse.  In  this  case,  most  of  the  zeros  that  were 
identified,  were  identified  incorrectly.  This  is  because  in 
this  case,  a  gradient  search  on  40  parameters  was  performed. 
This  once  again  leads  to  the  conclusion  that  another  method 
of  zero  identification  is  needed. 

The  BLS  method  did  not  create  an  accurate  state  space 
model  for  either  the  truth  model,  or  the  PACOSS  DTA.  Even 
in  the  fairly  single  2  input/2  output,  3  mode  truth  model, 
it  was  only  able  to  identify  2  of  the  modes,  and  those  modes 
were  not  identified  very  accurately.  Additionally,  the  BLS 
method  demonstrated  in  the  truth  model  a  vulnerability  to 
measurement  errors  of  as  small  as  1%  RMS.  This,  and  the 
fact  that  it  was  uncUale  to  identify  any  of  the  PACOSS  DTA 
modes  accurately  suggests  that  this  method  may  not  be 
suitable  for  identification  of  systems  of  the  magnitude  of 
the  PACOSS. 

Recommendations 

The  primary  recommendation  is  to  continue  working  with 
the  ERA  and  other  time  domain  methods  to  develop  an  accurate 
system  model  of  the  PACOSS.  The  system  identification, 
literature  as  well  as  the  tiruth  model  in  this  thesis, 
demonstrate  that  the  ERA  method  is  very  capable  of 
accurately  identifying  the  system  modes.  However,  for 
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control  system  design,  one  needs  the  input  and  output 
matrices  as  well.  The  ERA  method  is  not  very  accurate  at 
predicting  these  matrices.  Additional  work  should  be  done 
to  identify  a  method  that  will  identify  the  A,  B,  and  C 
matrix  independently,  or  in  conjunction  with  ERA,  but  more 
accurately  than  PEM. 
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Appendix  A.  Truth  Model 

This  appendix  contains  the  time  domain  and  bode  plots 
for  both  the  ERA/PEM  and  the  BLS/PEM  generated  models  for 
the  truth  system.  On  all  graphs,  the  truth  system  is  the 
solid  line  while  the  model  is  the  dashed  line. 
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The  ERA/PEM  methcxi  was  able  to  generate  a  state  space 
model  for  the  PACOSS  DTA.  Ihe  model  is  in  discrete  block 
diagonal  form,  and  contains  real  plant,  input,  and  output 
matrices.  The  plant  matrix  is  expressed  in  real,  block 
diagonal  form. 
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-0.0151 

0.0123 

-0.0039 

0.0164 

0.0173 
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-0.0729 

-0.0351 

0.3088 

0.1679 

0.1361 

0.0579 

0.68% 

2.1689 

-0.0691 

-0.0585 

-0.2660 

-0.0256 

-0.1183 

-0.0431 

0.4115 

-0.7177 

-0.1235 

0.0305 

0.0610 

-0.1052 

0.0442 

0.0643 

0.0113 

0.0941 

-0.0096 

0.0026 

-0.0119 

0.0144 

-0.0103 

0.0008 

0.3421 

0.1181 

-0.01454 

-0.0289 

0.0715 

0.0839 

-0.0743 

0.1312 

0.0037 

-1.1872 

0.0391 

0.0134 

-0.0212 

-0.0662 

0.0599 

0.0081 

-0.0215 

0.5127 

0.0607 

0.0369 

-0.0472 

-0.0275 

-0.2812 

-0.2208 

0.4714 

0.1425 

0.1203 

0.0372 

-0.1486 

-0.1580 

-0.0031 

-0.0186 

-0.4291 

0.0504 

-0.0611 

-0.0241 

-0.0363 

-0.0613 

0.0364 

0.0363 

-0.3400 

-0.0427 

-0.3156 

-0.0214 

-0.0402 

-0.0914 

0.0333 

0.0257 

-0.0805 

-0.0072 

-0.0195 

-0.3803 

-0.7671 

0.0800 

0.0603 

■0.1987 

-0.1103 

0.4800 

0.3328 

-0.0176 

0.0539 

-0.2555 

-0.4385 

■0.2449 

0.0299 

-0.1354 

0.0242 

-0.0050 

-0.0138 

0.0193 

-0.0319 

-0.0217 

-0.1942 

-0.0144 

0.0114 

-0.0022 

-0.0082 

0.0112 

-0.0144 

-0.0083 

-0.0817 

-0.0103 

-0.0880 

-0.0679 

0.0441 

-0.0235 

0.0497 

0.0134 

0.0152 

-0.0162 

0.0110 

0.0546 

-0.0137 

0.0090 

-0.0248 

-0.0106 

-0.0142 

0.0'T18 

0.6202 

-0.1401 

0.4226 

-0.3691 

-1.0326 

0.3141 

0.1365 

-0.2453 

0.2901 

-0.0494 

0.4162 

-0.3517 

-0.6825 

0.2415 

0.0803 

-1.9738 

■0.4932 

0.2255 

-♦.4244 

0.5451 

-1.4755 

0.0061 

■0.3389 

0.3047 

-0.6307 

0.07J8 

-0.5193 

-0.4821 

0.5788 

0.3518 

-0.5922 

0.0118 
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Appendix  C.  Bode  Plots 


This  appendix  contains  the  bode  plot  con^arissons  of 
the  actual  PACOSS  system  and  the  ERA/PEM  developed  model. 
All  plots  are  from  1  to  10  Hz,  and  are  frequency  response 
functions  of  velocity  output  to  force  input.  In  all  plots, 
the  dashed  line  represents  the  ERA/PEM  model  derived 
transfer  function,  while  the  solid  line  represents  the 
actual  PACOSS  DTA  derived  frequency  response  function. 
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This  appendix  contains  the  6  MATLAB  script  files  used 
in  this  thesis.  The  first  file  was  used  to  determine  the 


accuracy  of  the  ERA  model  by  computing  system  modes  euid 
comparing  the  results  to  those  of  the  truth  system.  The 
second  file  performed  the  same  euialysis  using  the  BLS  method 
to  determine  the  system  modes.  The  third  file  used  the  ERA 
method  to  conqpute  the  PACOSS  system  modes.  The  fourth  file 
used  the  BLS  method  to  ccm^jute  the  PACOSS  system  modes.  The 
fifth  computed  the  output  matrix  given  the  plant  auid  input 
matrices  from  either  the  ERA  method  or  the  BLS  method.  The 
final  file  confuted  and  con^ared  the  transfer  fuuictions  frcxn 
the  PACOSS  and  the  model  developed  by  either  the  ERA  or  the 
BLS  method. 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  erai 

%  This  file  generates  the  truth  information 

%  for  a  2  input  -  2  output  system  with  known 
%  A,  B,  and  C  matrices,  and  then 

%  identifies  that  system  using  the 

%  ERA  technique . 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  generate  input 
rand ( ' normal ' ) ; 

T  =  0.07; 
t  =  [0:T:20] ; 

A  =  1; 

%  generate  mesurement  errors 
E  =  0.15; 

for  1  =  l;length(t), 
e(l,l)  =  E*rcmd; 
e(l,2)  =  E/2*rand; 
end; 

%  Define  actual  system  TF 
%  Truth  Model 

%  Mass  coefficients 
ml»0. 61743;  m2=0. 19957;  m3=1.6582; 

%  Danqping  coefficients 

Cl=0. 031633;  c2=0. 32738;  C3=0. 32918;  c4=0. 019153; 
%  Stiffness  coefficients 

kl=52.058;  k2=0. 89183;  k3=19.756;  k4=178.08; 

%  Inverse  mass  matrix 
minv=diag( [l/ml  l/m2  l/m3] ) ; 

%  Stiffness  matrix 

kk=(kl+k2  -k2  0;-k2  k3+k2  -k3;0  -k3  k3+k4l ; 

%  Damping  matrix 

cc=[cl+c2  -C2  0;-c2  c3+c2  -C3;0  -C3  c3+C4] ; 

%  The  A,B,C,D  matrices  for  the  truth  system 
ag= (zeros (3, 3)  eye(3) ; -minv*kk  -minv*cc] ; 
bg*  (0  0;0  0;0  0;minv*[0  0;1  0;0  111; 
cg=[0  1  0  0  0  0;0  0  1  0  0  01  ; 
dg=zero8 (2, 2) ; 

%  Convert  to  a  TP 

[nuntgi,dengl]  «  ss2tf  (ag,bg,cg.dg,  l)  ; 

(nvung2,  deng2]  =  882tf  (ag,bg,  cg,dg,  2)  ; 
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[wn,  z]  =dainp(ag)  ; 

%  generate  bode  plot  of  truth  system 
w=logspace(-l,2,l50) ; 

[magi,phasl] =bode(ag,bg,cg,dg, l,w) ; 

[mag2,phas2] =bode{ag,bg, cg,dg,2,w) ; 

%  Now  change  to  discrete  state  space  form 
[ad,bd,cd,dd]  =  c2dm(ag,bg, cg,dg,T, ' zoh' ) ; 

%  generate  truth  output  (in^ulse  response  for  ERA) 
[yi,x]  =  impulse (ag,bg, eg, dg, 1, t) ; 

[y2,x]  =  impulse(ag,bg,cg,dg,2,t); 
yl (l;length{t) -1, : )  =  yl(2:length(t),:); 
y2 (l:length(t) -1, : )  =  yl(2:length(t) , : ) ; 

%  add  measurement  error  to  truth  output  signal 
for  1  =  l:length(t), 

yl(l,l)  =  yl(l,l)  +  e(l,l); 

yl{l,2)  =  yl{l,2)  +  e(l,2); 

y2{l,l)  =  y2(l,l)  +  e(l,l); 

y2(l,2)  =  y2(l,2)  +  e(l,2); 

end; 

y  =  yl+y2; 
yl=yl'; 
y2  =  y2 ' ; 

%  now  begin  system  identification 
ncols  =20; 
nrows  =60; 
inputs  =  2; 

[Y]  =  weave (yl,y2) ; 

[fdk, zmk, shapesk,partfak, EMACk, sv, at,bt, ct]  = 
erat (Y, l/T, ncols, nrows, inputs) ; 
cut  =  10; 
dt  =  zeros (2, 2); 

%  evaluate  model 
wdm=2  *pi *  f dk ; 
for  1  =  1 : length ( fdk) , 

wnm(l) =wdm(l) /sqrt (l-zmk(l) /100*zmk(l) /lOO) ; 
end; 

%  generate  random  input  for  PEM 
for  1  =  1 : length (t), 
u(l,l)  =  A*rand; 
u(l,2)  =  A*rand; 
end; 

%  generate  truth  output  from  rcuidom  input 
%  for  PEM 

[y,x]  *  lsim(ag,bg,cg,dg,u, t) ; 
y  »  y  +  e; 


74 


[junk, ninput] =size (bt) ; 

[noutput, junk]=size(ct) ; 
din=  zeros  (noutput,  ninput)  ; 

%  convert  bm  from  continuous  to  discrete 
[atemp,Bm]  =  c2d(at,bt,T) ; 

%  now  use  PEM  to  improve  C  one  row  at 
%  a  time  to  determine  zeros 
[rowa,cola]  =  size (Am); 

[rowb,colb]  =  size{Bm); 
ai  =  at; 
bi  =  bt; 

ci  =  nan*ones(l,rowa) ; 

di  =  zeros (1, colb) ; 

ki  =  zeros (rowa, 1) ; 

xOi  =  zeros (rowa, 1) ; 

msi  =  modstruc (ai,bi, ci,di,ki,xOi) ; 

parval  =  ct (l, : ) ; 

parva2  =  ct (2, : ) ; 

lambdi  =  []  ; 

thil  =  ms2th (msi, 'd' , parval, lambdi, T) ; 
thi2  =  ms2th(msi, 'd' ,parva2, lambdi, T) ; 
index  =  [1 : length (parval) ] ; 

%  perform  system  ID 

thl  =  pem([y(;,l)  u] , thil, index, -1, le-10, -1, -l, 
th2  =  pem([y(:,2)  u] , thi2, index, -1, le-10, -1, -1, 
%  convert  theta  to  SS  format 
[Am,Bm, cml,dm,kml,x01]  =  th2ss(thl); 

[Am,Bm,  cm2,dm,kmi,x01]  =  th2ss(th2); 

%  combine  c  matrices  to  get  MIMO  C  matrix 
Cm  =  [cmi;cm2l ; 

Dm  =  dt; 

for  1  =  l:length(t), 
u(l,l)  =  A*rand; 
u(l,2)  =  A*rand; 
end; 

[ya,x]  =  lsim(ag,bg,cg,dg,u, t) ; 

[ym,x]  =  dlsim(Am,  Bm,  Cm,Dm,  u) ; 

%  generate  bode  plot  of  model 

[magmi ,  phasml ]  =dbode  (Am,  Bm,  cml ,  dml ,  T,  1 ,  w)  ; 

[magin2 ,  phasm2  ]  =dbode  ( Am,  Bm,  cml ,  dml ,  T,  2 ,  w) ; 
plot(t,ya( : ,1) ,t,yml(; ,1) ) ; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  tlst 

%  This  file  generates  the  truth  information 
%  for  a  2  input  -  2  output  system  and  then 
%  identifies  that  system  using  a  backward 

%  batch  least  squares  technique.  The  results 

%  are  then  run  through  PEM  to  determine  the  C 

%  matrix 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  generate  input 
rand ( ' normal ' ) ; 

T  =  0.08; 
t  =  [0:T:20] ; 

A  =  1; 

%  generate  measurment  error 
S  =  0.01; 

for  1  =  l;length(t) -1, 
e(l,l)  =  E*rand; 
e(l,2)  =  E/2*rand; 
end; 

%  Define  actual  system  TP 
%  Truth  Model 

%  Mass  coefficients 
ml=0. 61743;  m2=0. 19957;  m3=1.6582; 

%  Damping  coefficients 

cl=0. 031633;  c2=0. 32738;  c3=0. 32918;  C4=0. 019153; 
%  Stiffness  coefficients 

kl=52.058;  k2=0. 89183;  k3=19.756;  k4=178.08; 

%  Inverse  mass  matrix 
minv=diag( [l/ml  l/m2  l/m3] ) ; 

%  Stiffness  matrix 

kk=[kl+k2  -k2  0;-k2  k3+k2  -k3;0  -k3  k3+k4] ; 

%  Daitping  matrix 

CC=  [Cl+C2  -C2  0;-c2  C3+c2  -C3;0  -C3  C3+C4] ; 

%  The  A,B,C,D  matrices  for  the  truth  system 
ag= [zeros(3,3)  eye{3) ; -minv*kk  -minv*cc] ; 
bg=[0  0;0  0;0  0;minv*[0  0;1  0;0  1]  1  ; 
cg=[0  1  0  0  0  0;0  0  1  0  0  0]; 
dg= zeros (2,2) ; 

%  Convert  to  a  TP 

[nijmgl,dengi]  =  ss2tf  (ag,bg,  cg,dg,  l) ; 
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[nunig2,deng21  =  ss2tf  {ag,bg, cg,dg,2)  ; 

[wn,  z]  =damp{ag)  ; 

%  generate  bode  plot  of  truth  system 
w=logspace{-l,2, 150) ; 

[magi,phasl] =bode {ag,bg, cg,dg, 1, w) ; 

[mag2,phas2] =bode (ag, bg, eg, dg, 2, w) ; 

%  Now  change  to  discrete  state  space  form 
[ad, bd, cd, dd]  =  c2dm(ag,bg, cg,dg,T, ' zoh' ) ; 

%  generate  random  inputs 
for  1  =  l:length(t), 
u{l,l)  =  A*rand; 
u{l,2)  =  A*rand; 
e(l,l)  =  E*rand; 
e(l,2)  =  E/2*rand; 
end; 

%  generate  truth  system  response  to 
%  random  input 
[y,x]  =  lsim(ag,bg,cg,dg,u, t) ; 

%  add  measurement  error  into  the  system 
y(l:l,:)  =y{l:l,:)  +e(l:l,:); 

%  now  to  identify  the  system 
p  =  6;  %  model  order 

[fd, zm,  zpoles, shapes, partfac]  =  mimo(u,y,p, l/T) 

%  find  the  undamped  natural  f regencies 

wdm=2*pi*fd; 

for  1  =  1 : length ( fd) , 

wnm(l) =wdm(l) /sqrt (l-zm(l) /l00*zm(l) /lOO) ; 
end; 

wnm  =  wnm'  ; 

%  create  the  state  space  model 

n= length (fd) ; 

nc=l; 

for  jh=l:n, 

wdl=2*pi*fd( jh) ; 

wnl=wdl/sqrt {l-zm( jh) /l00*zm( jh) /lOO) ; 
zKnc)  =exp(  (-zm(jh)  /100*wnl+j*wdl)  *T)  ; 
zl  (nc+1)  =conj  (zKnc) )  ; 
bm(nc, : ) =partfac{ jh, : ) ; 
bm(nc+l, : ) =conj (partfac (jh,  : ) )  ; 
cm( : , nc) =shapes ( : , jh) ; 
cm( : ,nc+l) =conj (shapes ( : , jh) ) ; 
nc=nc+2 ; 
end; 

am=diag(zl)  ; 

[junk,ninput] =size(bm) ; 
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[noiitput,  junk]  =83  26  (cm)  ; 
dm= zeros (noutput, ninput) ; 

%  convei't  bm  from  continuous  to  discrete 
[atemp,Bm]  =  c2d{am,bm,T) ; 

%  convert  am,  Bm,  and  cm  to  real  element  matrices 
T1  ..  [1  l;j  -j]; 

T2  -  Tl; 

[rl,cl]  =  size (am); 
for  k  =  4:2:rl, 
tr,c]  =  size{T2); 

T2  =  [T2  zeros {r,  2)  ;  zeros (2,  c)  Tl] ; 
end; 

Tinv  =  inv{T2) ; 

Am  =  T2*am*Tinv; 

Bm  =  T2*Bm; 

Cm  =  cm* Tinv, • 

Dm  =  dm; 

%  now  use  PEM  to  improve  C 
[rowa,cola]  =  size (Am)  ; 

[rowb,colb]  =  size(Bm); 

[atemp,Bt]  =  c2d (Am, Bm, T) ; 
ai  =  Am; 
bi  =  Bt; 

ci  =  nan*ones (1, rowa) ; 

di  =  zeros (1, colb) ; 

ki  =  zeros (rowa, 1) ; 

xOi  =  zeros (rowa, 1) ; 

msi  =  modstiruc(ai,bi,ci,di,ki,xOi) ; 

pairvai  =  Cm  ( i ,  : )  ; 

pairva2  =  Cm  ( 2 ,  : )  ; 

lambdi  =  []  ; 

thil  =  ms2th (msi, 'd' ,parval, lambdi, T) ; 
thi2  =  ms2th (ms i, 'd* ,parva2, lambdi, T) ; 
index  =  [1 : length (parval) ] ; 

%  perform  system  ID 

thl  =  pem([y(:,l)  u] , thil, index, -1, le-10, -1, -1,T) ; 
th2  =  pem([y(;,2)  u] , thi2, index, -l, le-10, -l, -1,T) ; 

%  convert  theta  to  SS  format 
[Am,Bm,  cml,dm,kml,x01]  =  th2ss(thl); 

[Am,Bm,  cm2,dm,  kml,x011  =  th2ss{th2); 

%  combine  c  matrices  to  get  MIMO  Cm  matrix 
Cm  =  [cml;cm2]  ; 

%  generate  bode  plot  of  system  model 
w=logspace (-1,2, 150) ; 

[magmi ,  phasml]  =dbode  (Am,  Bt ,  Cm,  Dm,  T,  l ,  w) ; 
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[niagin2,pha8in2]sdbode(Am,Bt,Cm,Dm,T,2,w)  ; 
%  conpare  results 
for  1  *  l:length(t), 
u(l,l)  ■  A*rand; 
u(l,2)  «  A*rand; 
end; 

[ya.x]  «  lsiin(ag,bg,cg,dg,u,t)  ;♦ 

[ym.x]  *  dlslm(Ain,Bm,Qn,Om,u); 
plot  (t,ya{  :,l),t,yin(:,l)) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  teral 

%  This  file  computes  the  system  model  for 
%  the  PACOSS  DTA  using  the  ERA  technique 
%  It  uses  swept  sine  data  to  generate  acceleration 
%  transfer  functions,  then  integrates  those  values 
%  to  obtain  velocity  TFs,  and  then  takes 
%  the  inverse  fourier  transform  to  obtain 

%  the  impulse  response.  The  impulse  response 

%  is  then  inserted  into  ERA  to  obtain  the  system 
%  modes . 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  load  the  response  data  and  create  the  y  vectors 
%  input  #1 


load 

aiml 

load 

aim2 

load 

aim3 

load 

aim4 

load 

aims 

load 

alm6 

load 

alm9 

load 

almlO 

load 

a2ml 

load 

a2m2 

load 

a2m3 

load 

a2m4 

load 

a2m5 

load 

a2m6 

load 

a2m9 

load 

a2ml0 

load 

a3ml 

load 

a3m2 

load 

a3m3 

load 

a3m4 

load 

a3m5 

load 

a3m6 

load 

a3m9 

load 

a3ml0 

%  input  #2 


%  input  #3 


%  input  #4 
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load  a4mi 
load  a4in2 
load  a4m3 
load  a4m4 
load  a4m5 
load  a4m6 
load  a4m9 
load  a4inl0 

%  input  #5 

load  a5ml 
load  a5ni2 
load  a5m3 
load  a5m4 
load  a5m5 
load  a5m6 
load  a5m9 
load  a5ml0 

%  input  #6 

load  a6ml 
load  a6m2 
load  a6m3 
load  a6m4 
load  a6m5 
load  a6m6 
load  a6m9 
load  a6ml0 

%  input  #9 

load  a9ml 
load  a9in2 
load  a9ni3 
load  a9ra4 
load  a9m5 
load  a9m6 
load  a9m9 
load  a9ml0 

%  input  #10 

load  alOml 
load  al0in2 
load  ai0in3 
load  al0m4 
load  alOmS 
load  al0m6 
load  al0m9 
load  alOmlO 

%  define  the  frequency  vector 
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freq  =  0:0.025:10; 

%  now  the  transfer  function  data  must  be  recoiranbined  into 
%  con^lex  numbers . 
for  1  =  1:401, 

%  input  #1 

aii(l)  =  almi(l)  +  j*aimi (1+801) ; 
al2(l)  =  alm2{l)  +  j*alm2 (1+801) ; 
al3(l)  =  alm3(l)  +  j*alm3  (1+801) ; 
ai4(l)  =  aim4(l)  +  j*alm4 (1+801) ; 
ai5(l)  =  aimsd)  +  j*alm5  (1+801)  ; 
al6(l)  =  aimed)  +  j*alm6 (1+801)  ; 
al9(l)  =  alm9(l)  +  j*alm9  (1+801) ; 
allOd)  =  almlOd)  +  j*alml0  (1+801)  ; 

%  input  #2 

a2i(l)  =  a2mi(l)  +  j*a2ml (1+801) ; 
a22(l)  =  a2in2(l)  +  j*a2m2  (1+801) ; 
a23(l)  =  a2m3(l)  +  j*a2m3 (1+801) ; 
a24(l)  =  a2m4(l)  +  j*a2m4 (1+801) ; 
a25(l)  =  a2m5(l)  +  j*a2m5 (1+801) ; 
a26(l)  =  a2m6(l)  +  j*a2m6 (1+801) ; 
a29(l)  =  a2m9d)  +  j*a2m9  (1+801)  ; 
a2i0(l)  =  a2ml0(l)  +  j*a2ml0 (1+801) ; 

%  input  #3 

a31(l)  =  a3mi(l)  +  j*a3mld+801)  ; 
a32d)  =  a3m2d)  +  j*a3m2  (1+801)  ; 
a33(l)  =  a3m3(l)  +  j*a3m3  (1+801) ; 
a34(l)  =  a3m4d)  +  j*a3m4  (1+801) ; 
a35(l)  =  a3m5(l)  +  j*a3m5  (1+801) ; 
a36d)  =  a3m6(l)  +  j*a3m6 (1+801)  ; 
a39(l)  =  a3m9(l)  +  j*a3m9  (1+801) ; 
a310(l)  =  a3ml0d)  +  j*a3ml0  (1+801)  ; 

%  input  #4 

a4l(l)  =  a4mi(l)  +  j *a4ml (1+801) ; 
a42(l)  =  a4in2(l)  +  j*a4m2  (1+801)  ; 
a43(l)  =  a4m3d)  +  j*a4m3  (1+801)  ; 
a44(l)  =  a4m4(l)  +  j*a4m4  (1+801) ; 
a45(l)  =  a4in5d)  +  j*a4m5 (1+801) ; 
a46(l)  =  a4m6(l)  +  j*a4m6 (1+801) ; 
a49(l)  =  a4m9(l)  +  j*a4m9 (1+801) ; 
a4l0(l)  =  a4ml0(l)  +  j*a4ml0 (1+801) ; 

%  input  #5 

aSld)  =  a5ml(l)  +  j*a5ml  (1+801)  ; 
a52(l)  *  a5m2(l)  +  j*a5m2 (1+801) ; 
a53(l)  =  a5m3(l)  +  j*a5m3 (1+801) ; 
a54(l)  =  a5m4d)  +  j*a5m4  (1+801); 
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a55(l)  =  aSmSd)  +  j*a5m5 (1+801)  ; 
a56(l)  =  aSmSd)  +  j*a5m6 (1+801) ; 
a59(l)  =  a5m9(l)  +  j*a5m9 (1+801) ; 
a510(l)  =  aSmlOd)  +  j*a5ml0  (1+801)  ; 

%  input  #6 

a61(l)  =  aemld)  +  j*a6ml (1+801)  ; 
a62(l)  =  a6m2(l)  +  j*a6m2 (1+801) ; 
a63(l)  =  a6m3(l)  +  j*a6m3  (1+801) ; 
a64(l)  =  a6m4(l)  +  j*a6m4 (1+801) ; 
a65(l)  =  aemSd)  +  j*a€m5 (1+801)  ; 
a66(l)  =  a6m6(l)  +  j *a6in6 (1+801)  ; 
a69(l)  =  a6iti9(l)  +  j*a6m9  (1+801)  ; 
aGlOd)  =  aemiOd)  +  j*a6ml0  (1+801)  ; 

%  input  #9 

a91(l)  =  a9inl(l)  +  j*a9ml (1+801)  ; 
a92(l)  =  a9tn2(l)  +  j*a9m2  (1+801)  ; 
a93(l)  =  a9m3(l)  +  j*a9m3 (1+801) ; 
a94(l)  =  a9m4(l)  +  j*a9m4 (1+801) ; 
a95(l)  =  a9m5(l)  +  j*a9m5  (1+801) ; 
a96(l)  =  a9m6(l)  +  j*a9m6  (1+801) ; 
a99(l)  =  a9m9(l)  +  j*a9in9  (1+801)  ; 
a9l0(l)  =  a9ml0(l)  +  j*a9ml0 (1+801) ; 

%  input  #10 

alOld)  =  alOmld)  +  j*al0ml (1+801) ; 
ai02(l)  ®  al0m2(l)  +  j*alOin2 (1+801) ; 
al03(l)  =  al0m3(p  +  j*al0m3  (1+801) ; 
al04(l)  =  alOm4(i)  +  j*al0m4 (1+801) ; 
ai05(l)  =  alOmSd)  +  j*al0in5  (1+801) ; 
al06(l)  =  al0m6(l)  +  j*al0in6 (1+801)  ; 
ai09(l)  =  alOm9(l)  +  j*al0m9 (1+801) ; 
alOlOd)  =  alOmlOd)  +  j*alOmlO (1+801)  ; 
end; 

%  integrate  to  get  velocity 
%  set  0  Hz  value  to  zero  (no  rigid  body  modes) 
%  input  #1 

vll(l)  =  0;vl2(l)  =  0;vl3(l)  =  0; 

vl4(l)  =  0;vl5(l)  =  0;vl6(l)  =  0; 

V19(1)  =  0;vll0(l)  =  0; 

%  input  #2 

V21(l)  =  0;v22(l)  =  0;v23(l)  =  0; 

V24(l)  =  0;v25(l)  =  0;v26(l)  =  0; 

v29(l)  =  0;v210(l)  =  0; 

%  input  #3 

v31(l)  =  0;v32(l)  =  0;v33(l)  «  0; 

v34(l)  *  0;v35(l)  =  0;v36(l)  »  0; 
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v39(l)  «  0;v310(l)  »  0; 

%  input  #4 

V41(l)  »  0;v42{l)  «  0;v43(l)  =  0; 

V44(l)  a  0;V45(1)  *  0;v46(l)  =  0; 

v49(l)  =  0;v410(l)  =  0; 

%  input  #5 

v51(l)  =  0;v52(l)  =  0;vl3(l)  =  0; 

v54(l)  s  0;V55(1)  *  0;vl6(l)  =  0; 

V59(l)  =  0;v510(l)  =  0; 

%  input  #6 

v61(l)  =  0;v62(l)  *  0;vl3(l)  -  0; 

v64(l)  =  0;v65(l)  *  0;vl6(l)  »  0; 

v69(l)  =  0;v610(l)  -  0; 

%  input  #9 

v91(l)  =  0;v92(l)  *  0;vl3(l)  *  0; 

v94(l)  =  0;V95(1)  »  0;vl6(l)  «  0; 

v99(l)  =  0;v910(l)  »  0; 
t  input  #10 

vlOl(l)  «  0;V102(1)  -  0;vl3(l)  -  0; 
vl04(l)  «  O.'VlOSd)  -  0;vl6(l)  -  0; 
vl09(l)  -  0;vl010(l)  «  0; 

%  divide  by  jw  to  integrate 
for  1  «  3:401, 

8  «  j*2*pi*freq(l) ; 

%  input  #1 
V11{1,1)  .  all(l)/8; 
vi2(l,l)  -  al2(l)/8; 
vi3(l,l)  -  al3(l)/8; 

VM(1,1)  -  al4(l)/8; 
vi5(l,l)  -  al5(l)/8; 
vi6(l, 1)  -  al6(l) /8; 
vi9(l, 1)  -  ai9(l) /8; 
vllOd.l)  .  all0(l)/8; 

%  input  #2 
v21{l,l)  -  a21(l)/8; 
v22(l,l)  .  a22(l) /8; 
v23 (1 . 1)  -a23(l)/B; 

V24<1. 1)  -  824(1) /8; 

v2S(l. 1)  -  825(1) /8, 
v24 (1.1)  -  824(1) /8, 
v29(l, 1)  -  829(1) /8, 

V210(1, 1)  .  8210(1) /8. 

%  input  #3 
V31  (1,  1)  .  831  (1)  /B, 
v32 (1.1)  -  832 ( 1 ) /a. 
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v33(l,l)  =  a33(l)/s; 
v34(l,l)  =  a34(l)/s; 
v35(l,l)  =  a35(l)/s; 
v36{l,l)  =  a36(l)/s; 
v39(l,l)  =  a39{l)/s; 
v310(l,l)  =  a310(l)/8; 

%  input  #4 
v41(l,l)  =  a41{l)/s; 
V42(l,l)  =  a42(l)/S; 
V43(l,l)  =  a43(l)/8; 
V44(l,l)  =  a44(l)/8; 
V45(l,l)  *=  a45{l)/8; 
V46(l,l)  =  a46{l)/8; 
v49(l,l)  =  a49(l)/8; 
V410(l,l)  «  a410(l)/8; 

%  input  #5 
v51(l,l)  »  a51(l)/8; 
V52(l,l)  -  a52(l)/8; 
v53(l,l)  «  a53(l)/8; 
V54(l,l)  *  a54(l)/8; 
V55(l,l)  »  a5S(l)/8; 
V56(l,l)  -  a56(l)/8; 
V59(l,l)  •  a59(l)/8; 
V510(l,l)  -  a510(l)/s; 

%  input  #6 
V61(l,l)  -  a61(l)/8; 
V62 (1, 1)  *  a62 (1) /s; 
v63 (1, 1)  «  a63 (1) /s; 
V64 (1, 1)  «  a64 (1) /s; 
V65(l, 1)  -  a65(l) /8; 
V66(l, 1)  -  a66(l) /b; 

v69 (1, 1)  -  a69 (1) In, 
V610(l,l)  -  a610(l)/8; 

%  input  19 
V91 (1, 1)  -  a91 (1) /8; 
V92(l. 1)  -  a92(l) /8. 
V93  ( 1,1)  -  a93  <  1  )  /8, 

v94 (1,1)  -  a94 ( 1 )  /t , 
V9S(1, 1)  -  a9Sll) /a. 
V94 (1,1)  •  a94 ( 1 ) /b, 
v99 (1,1)  -  a99 ( 1 ) /B, 
v^io ( 1 .  1 )  •  a9io ( 1 ) /B, 

%  input  910 
vlOl (1,1)  -  alOl ( 1 ) /B, 
Vloa (1,1)  -  ai02 ( 1) /B, 


vl03(l,l)  *  al03{l)/8; 
vl04(l,l)  =  al04(l)/8; 
vl05(l,l)  =  al05(l)/8; 

V106(1,1)  *  al06(l)/S; 
vl09(l,l)  =  al09(l)/8; 
vlOlOd.l)  =  al010(l)/8; 
end; 

%  find  0.025  Hz  value  (eet  to  small  number  so  it  does  not 
%  interfere  with  results) 

%  input  #1 

vll{2,l)  -  vll(3,l)  ;V12{2,1)  =  vl2(3.1)  ;vl3(2,l)  =vl6(3,l) 
vl4(2.1)  *  V14(3,1)  ;V15(2,1)  =  V15(3,1)  ;vl6(2,l)  =vl6(3,l) 
vl9(2,l)  =  V19(3,1) ;V110(2,1)  *vll0(3,l); 

%  input  #2 

v21(2,l)  *  V21(3,1) ;v22(2,l)  =  v22 (3. 1) ;v23 (2, 1)  =  v26(3,l) 
v24(2,l)  *  v24(3,l)  ;v25(2,l)  *  v25  (3. 1)  .*  V26  (2, 1)  =v26(3,l) 
v29(2,l)  »  v29(3,l) ;v210(2,l)  *v210(3,l); 

%  input  #3 

v31{2,l)  »  V31(3,l) ;v32(2,l)  =  v32 (3, 1) ;v33 (2 , 1)  =v36(3,l) 
v34(2,l)  «  v34(3,l)  ;V35(2,1)  =  v35  (3 , 1)  ;  v36  (2 , 1)  =v36(3,l) 
v39(2,l)  -  v39(3,l) ;v310(2,l)  =v310(3,l); 

%  input  #4 

V41(2.1)  «  V41(3,l) ;V42(2,1)  =  V42 (3, 1) ;v43 (2, 1)  »v46(3,l) 
V44(2,l)  •  v44(3,l) ;V45(2,1)  »  V45 (3, 1) ;V46 (2, 1)  =  V46{3,1) 
V49(2,l)  -  v49(3,l) ;v410(2,l)  *v410(3,l); 

%  input  #5 

v51(2,l)  »  v51 (3, 1) ;V52 (2, 1)  -  v52 (3, 1) ; v53 (2, 1)  =  V56{3,1) 

v54(2.1)  «  V54(3,l) ;V55(2.1)  -  v55 (3, 1) ; v56 (2. 1)  =  v56(3,l) 

v59(2,l)  *  v59(3, 1) ;v510(2.1)  *  v510(3,l); 

%  input  #6 

v61(2.1)  -  v61 (3, 1) ; v62 (2, 1)  *  v62 (3, 1) ; v63 (2, 1)  =  v66(3,l) 

v64(2,l)  -  v64  (3. 1)  ;v65(2, 1)  »  v€5  (3 , 1)  .*  v66  (2 , 1)  =  v66(3,l) 

v69(2,l)  -  v69(3,l) ;v610(2,l)  *v610(3,l); 

%  input  #9 

v91(2,l)  -  v91 (3, 1) ; v92 (2, 1)  -  v92 (3, 1) ; v93 (2, 1)  =  v96(3,l) 
V94(2.1)  -  v94 (3. 1) ;v95(2, 1)  -  v95 ( 3 , 1) ; v96 (2 , 1)  =v96(3,l) 
v99{2,l)  -  v99 (3, 1) ; v910 (2, 1)  »v910(3,l); 

%  input  810 

vl01(2,l)  .  vlOl (3. 1) ;vl02 (2. 1)  •  vl02 (3, 1) ;vl03 (2, 1)  » 

V10413. 1) , 

vl04(2.1)  »  V104  (  3 , 1 )  ;  vlOS  (2, 1)  •  Vl05  (  3 , 1 )  .' Vl06  (2 , 1)  * 

V106  (3.1)  . 

vl09(2.1)  -  V109 ( 3. 1) ; vlOlO (2, 1)  «  Vl0l0(3,l); 

%  due  to  poor  data,  reduce  influence  of  2nd  and  3rd  point 
for  1  -  23, 


vll(l,l)  =  O.Ol*vll{l,l) ;vl2(l, 1)  =  0.01*vl2{l,l) ;vl3{l,l) 
=  0.01*vl3 (1,1); 

Vl4(l,l)  =  0.01*V14{1,1) ;vl5(l,l)  =  0.01*V15{1, 1) ;vl6(l, 1) 
=  0.01*vl6(l,l) ; 

vl9(l,l)  =  0.01*vl9(l,l) ;vllO(l,l)  =  0 . 01*vll0 (1 , 1) ; 
v21(l,l)  =  0.01*v21{l,l) ;v22(l.l)  =  0 . 01*v22 (1, 1) ; v23 (1, 1) 
=  0.01*v23(l,l) ; 

v24{l.l)  :=  0.01*v24{l,l)  ;v25{l,l)  =  0 . 01*v25  (1 , 1)  ;  v26  (1 , 1) 
=  0 . 01*v26 (1, 1) ; 

V29{1,1)  =  0.01*v29(l,l) ;v210{l,l)  =  0 . 01*v210 (1, 1) ; 
v31(l,l)  =  0.01*v31(l,l) ;v32(l,l)  =  0 . 01*v32 (1, 1) ;v33 (1, 1) 
=  0 . 01*v33 (1,1)  ; 

v34(l,l)  =  0.01*v34(l,l) ;v35(l,l)  =  0 . 01*v35 (1, 1) ;v36 (1, 1) 
=  0.01*v36(l.l) ; 

v39(l,l)  =  0.01*v39{l.l) ;v310{l,l)  =  0 . 01*v310 (1, 1) ; 
v41{l,l)  =  0.01*v41(l,l) ;v42(l,l)  =  0.01*v42(l,l) ;v43(l,l) 
=  0.01*v43(l,l) ; 

V44{1,1)  =  0.01*V44(1,1) ;v45(l,l)  =  0 . 01*v45 (1, 1) ;v46 (1, 1) 
=  0.01*v46(l,l) ; 

v49(l,l)  =  0.01*v49 (1,1) ;v410{l,l)  =  0 . 01*v410 (1, 1> ; 
v51(l,l)  =  0.01*v51(l,l) ;v52{l,l)  =  0.01*v52(l,l) ;v53(l,l) 
=  0.01*v53(l.l)  ; 

v54{l,l)  =  0.01*v54(l,l) ;v55(l,l)  =  0 . 01*v55 (1, 1) ;  v56 (1, 1) 
=  0.01*v56(l,l)  ; 

v59(l,l)  =  0.01*v59(l,l) ;v510{l,l)  =  0 . 01*v510 (1, 1) ; 
v61(l,l)  =  0.01*v61(l,l) ;v62{l,l)  =  0.01*v62{l,l) ;v63(l,l) 
=  0.01*v63{l,l) ; 

v64(l,l)  =  0.01*v64(l,l) ;v65(l,l)  =  0 . 01*v65 (1, 1) ; v66 (1, 1) 
=  0 . 01*v66 (1, 1) ; 

v69(l,l)  =  0.01*v69(l,l) ;v610(l,l)  =  0 . 01*v610 (1, 1) ; 
v91{l,l)  =  0.01*v91(l,l) ;v92(l,l)  =  0 . 01*v92 (1, 1) ; v93 (1, 1) 
=  0.01*v93 (1, 1) ; 

v94(l,l)  =  0.01*v94(l,l) ;v95(l,l)  =  0 . 01*v95 (1, 1) ;v96 (1, 1) 
=  0.01*v96{l,l) ; 

v99{l,l)  =  0.01*v99(l,l) ;v910{l,l)  =  0 . 01*v910 (1, 1) ; 
vl01(l,l)  =  0.V)l*vl01{l,l)  ;vl02(l,l)  = 

0.01*vl02(l,l) ;vl03{l,l)  =  0 . 01*vl03 (1, 1) ; 

vl04(l,l)  =  0.01*vl04(l,l) ;vl05(l,l)  = 

0.01*vl05(l,l) ;vl06(l,l)  =  0 . 01*vl06 (1, 1) ; 

vl09(l,l)  =  0.01*vl09(l,l) ;vl010(l,l>  =  0 . 01*vl010 (1, 1) ; 
end; 

%  save  the  velocity  TF's 

save  tfvl.mat  vll  vl2  vl3  vl4  vl5  vl6  vl9  vllO 

save  tfv2.mat  v21  v22  v23  v24  v25  v26  v29  v210 

save  tfv3.mat  v31  v32  v33  v34  v35  v36  v39  v310 
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save  tfv4.mat  v41  v42  v43  v44  v45  v46  v49  v410 

save  tfvB.mat  v51  v52  v53  v54  v55  v56  v59  v510 

save  tfvS.mat  v61  v62  v63  v64  v65  v66  v69  v610 

save  tfv9.inat  v91  v92  v93  v94  v95  v96  v99  v910 

save  tfvlO.mat  vlOl  vl02  vl03  vl04  vl05  vl06  vl09  vlOlO 
%  tcJce  the  inverse  fourier  transform  to  get  the 
%  impulse  response  (2048  point  of  data  were  gathered) 
%  input  #1 

yll  =  2048*real(ifft(vll,2048) ) ; 
yl2  =  2048*real(ifft(vl2,2048) ) ; 
yl3  =  2048*real{ifft(vl3,2048) ) ; 
yl4  =  2048*real(ifft{vl4,2048) ) ; 
yl5  =  2048*real(ifft(vl5,2048) ) ; 
yl6  =  2048*real(ifft(vl6,2048) ) ; 
yl9  =  2048*real (ifft (vl9,2048) ) ; 
yllO  =  2048*real (ifft (vllO, 2048) ) ; 

%  input  #2 

y21  =  2048*real (ifft (v21, 2048) ) ; 
y22  =  2048*real (ifft (v22, 2048) ) ; 
y23  =  2048*real(ifft(v23,2048) )  ; 
y24  =  2048*real (ifft (v24, 2048) ) ; 
y25  =  2048*real (ifft (v25, 2048) )  ; 
y26  =  2048*real ( ifft (v26, 2048) ) ; 
y29  =  2048*real (ifft (v29, 2048) )  ; 
y210  =  2048*real (ifft (v210, 2048) ) ; 

%  input  #3 

y31  =  2048*real (ifft (v31, 2048) ) ; 
y32  =  2048*real (ifft (v32, 2048) )  ; 
y33  =  2048*real (ifft (v33, 2048) ) ; 
y34  =  2048*real (ifft (v34, 2048) ) ; 
y35  =  2048*real (ifft (v35, 2048) ) ; 
y36  =  2048*real (ifft (v36, 2048) ) ; 
y39  =  2048*real (ifft (v39, 2048) ) ; 
y3l0  =  2048*real (ifft (v310, 2048) ) ; 

%  input  #4 

y41  =  2048*real (ifft (v41, 2048) ) ; 
y42  =  2048*real (ifft (v42, 2048) ) ; 
y43  *  2048*real (ifft (v43, 2048) ) ; 
y44  =  2048*real (ifft (v44, 2048) ) ; 
y45  =  2048*real (ifft (v45, 2048) ) ; 
y46  =  2048*real (ifft (v46, 2048) ) ; 
y49  =  2048*real (ifft (v49, 2048) ) ; 
y410  =  2048*real(ifft(v410,2048) ) ; 

%  input  #5 

ySl  =  2048*real (ifft (v51, 2048)); 
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y52  »  2048*real ( Iff t(v52, 2048) ) ; 
y53  »  2048*real(lfft(v53.2048) ) ; 
y54  =  2048*real(ifft(v54,2048) ) ; 
y55  =  2048*real(ifft(v55,2048) ) ; 
y56  =  2048*real(ifft(v56,2048) )  ; 
y59  =  2048*real{ifft(v59,2048) )  ; 
ySlO  =  2048*real (ifft (v510,2048) ) ; 

%  input  #6 

y61  =  2048*real (Ifft (v61, 2048) ) ; 
y62  =  2048*real (ifft (v62, 2048) ) ; 
y63  »  2048*real (ifft (v63, 2048) ) ; 
y64  =  2048*real (ifft (v64, 2048) )  ; 
y65  =  2048*real (ifft (v65, 2048) ) ; 
y66  =  2048*real (ifft (v66, 2048) ) ; 
y69  =  2048*real (ifft (v69, 2048) ) ; 
ySlO  =  2048*real (ifft (v610, 2048) ) ; 

%  input  #9 

y91  =  2048*real (ifft (v9l, 2048) ) ; 
y92  =  2048*real (ifft (v92, 2048) ) ; 
y93  =  2048*real (ifft (v93, 2048) ) ; 
y94  =  2048*real (ifft (v94, 2048) ) ; 
y95  =  2048*real (ifft (v95, 2048) ) ; 
y96  »  2048*real (ifft (v96, 2048) ) ; 
y99  »  2048*real (ifft (v99, 2048) ) ; 
y910  =  2048*real (ifft (V910, 2048) ) ; 

%  input  #10 

ylOl  =  2048*real (ifft (vlOl, 2048) ) ; 
yl02  =  2048*real (ifft (vl02, 2048) ) ; 
yl03  =  2048*real (ifft (vl03, 2048) ) ; 
yl04  =  2048*real (ifft (vl04, 2048) ) ; 
yl05  =  2048*real (ifft (vl05, 2048) ) ; 
yl06  =  2048*real (if ft (vl06, 2048) ) ; 
yl09  =  2048*real (ifft (vl09, 2048) ) ; 
ylOlO  =  2048*real (ifft (vlOlO, 2048) ) ; 

%  remove  uneeded  vectors  to  clear  memory  space 

clear  vll  vl2  vl3  vl4  vl5  vl6  vl9  vllO 

clear  v21  v22  v23  v24  v25  v26  v29  v210 

clear  v31  v32  v33  v34  v35  v36  v39  v310 

clear  v41  v42  v43  v44  v45  v46  v49  v410 

clear  v51  v52  v53  v54  v55  v56  v59  v510 

clear  v61  v62  v63  v64  v65  v66  v69  v610 

clear  v9l  v92  v93  v94  v95  v96  v99  v910 

clear  vlOl  vl02  vl03  vl04  vl05  vl06  vl09  vlOlO 

%  conbine  the  incise  response  vectors  into  matrices 

zl  =  [yll  yl2  yl3  yl4  yl5  yl6  yl9  yllO] ; 
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z2  »  ty21  y22  y23  y24  y25  y26  y29  y210] ; 

z3  »  [y31  y32  y33  y34  y35  y36  y39  y310] ; 

z4  =  ty41  y42  y43  y44  y45  y46  y49  y410] ; 

z5  =  [y51  y52  y53  y54  y55  y56  y59  yBlO] ; 

z6  =  [y61  y62  y63  y64  y65  y66  y€9  ySlO] ; 

z9  =  [y91  y92  y93  y94  y95  y96  y99  y910] ; 

zlO  =  [ylOl  yl02  yl03  yl04  ylOS  yl06  yl09  ylOlO] ; 

%  remove  \ineeded  vectors  to  clear  memory  space 

clear  yll  yl2  yl3  yl4  yl5  yl6  yl9  yllO 

clear  y21  y22  y23  y24  y25  y26  y29  y210 

clear  y31  y32  y33  y34  y35  y36  y39  y310 

clear  y41  y42  y43  y44  y45  y46  y49  y410 

clear  y51  y52  y53  y54  y55  y56  y59  y510 

clear  y61  y62  y63  y64  y65  y66  y69  ySlO 

clear  y91  y92  y93  y94  y95  y96  y09  y910 

clear  ylOl  yl02  yl03  yl04  yl05  yl06  yl09  ylOlO 

%  change  from  row  to  column  form 

yl  =  zl • ; 

y2  =  z2  '  ; 

y3  =  z3 • ; 

y4  =  z4 ' ; 

yS  =  z5'; 

y6  =  z6 ' ; 

y9  =  z9'; 

ylO  =  zlO ' ; 

%  begin  system  identification 

ncols  =  90; 
nrows  =  400; 
inputs  =  8; 

%  set  sanple  rate 
fs  =  51.2; 

[Y]  =  weave(yl,y2,y3,y4,y5,y6,y9,yl0); 

[fd,zm, shapes, partfac,EMAC,sv]  » 
era (Y, fs, ncols, nrows, inputs) ; 
cut  =  0.0; 

%  eliminate  undesired  modes 

keep  =  input  Chow  many  modes  will  you  keep  '}; 
for  1  «  l:keep, 

k  a  input ( ' keep  mode  # :  ’ ) ; 

fdk(l,l)  =  fd(k); 

zmk(l, 1)  a  zm(k) ; 

shapesk ( : , 1)  a  shapes ( : , k) ; 

partfakd,:)  a  partfac(k,  :) ; 
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BMHOcd,!)  «  BBAC(k,l); 
end; 

%  evaluate  model 

wdm=2*pi*fdk; 

for  1  =  1 : length ( fdk) , 

wnm(l)  swdm(l)  /sqrt  (l-zmkd)  /I00*zmkd)  /lOO)  ; 
end; 

wnm  =  wnm/2/pl; 

%  Create  state  space  model  (A,b,C,  amd  D  matrices) 

[am,  bm,  cm,  dm]  » 

erastate(fdk,  zmk:,  shape sk, par t£ak,EMACk,  £s,  cut)  ; 

%  convert  hot  from  continuous  to  discrete 
[atenip,Bm]  =  c2d(am,to,  l/£s) ; 

%  convert  am,Bm,  and  cm  to  real  element  matrices  only 
T1  =  [1  l;j  -jl; 

T  *  Tl; 

[rl,ci]  «  slze(am); 
for  k  -  4:2:rl, 

[r,c]  *  slze(T); 

T  «  [T  zeros (r, 2) ; zeros (2, c)  Tl] ; 
end; 

Tlnv  «  Inv(T) ; 

A  »  T*am*Tlnv; 

B  «  T*Bm; 

C  »  cm*Tlnv; 

D  >  dm; 

save  tabc.mat  A  B  C  D 

save  tresul.mat  fd  zm  shapes  partfac  EMAC  sv 
save  tkept.mat  wnm  fdk  zmk  shapesk  partfaUc  EMACk 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  pera 

%  This  file  computes  the  system  model  for 
%  the  PACOSS  DTA  using  the  BRA  technique 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  load  the  bias  for  each  accelerometer 
load  noml2; 

blasd/  =  mean(Tlmel3) ; 
blas(2)  =  nieaui(Tlmel4) ; 
load  nom34; 

blasO)  =  me«ui(Tlmel3) ; 
blas(4)  =  mecm(Tlmel4)  ; 
load  ncmi56; 

blas(5)  =  mean(Tlmel3)  ; 
blas(6)  =  meeLn(Tlmel4)  ; 
load  nomSlO; 
blasO)  a  mean(Tlmel3)  ; 
blas(lO)  -  mean(Tlmel4) ; 

%  load  the  response  data  and  create  the  y  vectors 
%  Input  #1 
load  remll.mat; 
ul ( : , 1)  =  Timeil; 
yl ( : , 1)  =  Tiroei2; 
load  ranl2.mat; 
ul( : ,2)  =  Timeil; 
yl( : ,2)  =  Timei2; 
load  ranl3.mat; 
ul{:,3)  =  Timeil; 
yl  ( : ,  3 )  =  Tixnei2  ; 
load  rcuil4.mat; 
ul ( : , 4)  =  Timeil; 
yl( : ,4)  =  Timei2; 
load  reuiis.mat; 
ul{:,5)  =  Timeil; 
yl ( : , 5)  =  Timei2; 
load  ranl6.mat; 
ul ( : , 6 )  »  Timeil ; 
yl( : ,6)  *  Timei2; 
load  rani9.mat; 
ul( : ,7)  =  Timeil; 
yl( ; ,7)  *  Timei2; 
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%  input  #2 

load  ran21.mat; 

\i2  ( : ,  1)  =  Tlmell; 
y2  ( : ,  1)  =  Tiinei2; 
load  ran22.inat; 
u2  ( : ,2)  =  Tlmell; 
y2  ( : , 2 )  =  Tlmel2 ; 
load  rcui23.mat; 
u2 { : , 3 )  =  Tlmell ; 
y2  ( : , 3 )  =  Timei2 ; 
load  ran24.mat; 
u2(:,4)  =  Tlmell; 
y2(:,4)  =  Timei2; 
load  ran25.mat; 
u2(:,5)  =  Tlmell; 
y2 ( : , 5)  =  Timei2; 
load  ran26.mat; 
u2 ( : , 6)  =  Tlmell; 
y2{:,6)  =  Timei2; 
load  rauri29.mat; 
u2 ( : , 7)  =  Tlmell; 
y2 ( : , 7)  =  Timei2; 

%  input  #3 

load  ran31.mat; 
u3  ( : , 1)  =  Tlmell; 
y3  ( : , 1)  =  Timei2; 
load  ran32.mat; 
u3 ( : , 2)  =  Tlmell; 
y3(:,2)  «  Tlmel2; 
load  raui33.mat; 
u3  ( : , 3)  =  Tlmell; 
y3  ( : , 3)  =  Timei2; 
load  ran34.mat; 
u3(:,4)  =  Tlmell; 
y3  ( : ,4)  =  Timei2; 
load  ran35.mat; 
u3(;,5)  =  Tlmell; 
y3(:,5)  =  Timei2; 
load  ran36.mat; 
u3(: ,6)  =  Tlmell; 
y3( : ,6)  -  Timei2; 
load  ran39.mat; 
u3(: ,7)  -  Tlmell; 
y3 ( : , 7)  «  Timei2; 

%  Input  #4 
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load  reui41.inat; 
u4 ( : , 1)  =  Timeil; 
y4  ( : ,  1)  =  Tiinei2; 
load  rcui42.inat; 
u4(:,2)  =  Timeil; 
y4 ( : , 2)  =  Timei2; 
load  ran43.mat; 
u4(:,3)  =  Timeil; 
y4(:,3)  =  Timei2; 
load  ran44.mat; 
u4{:,4)  =  Timeil; 
y4(:,4)  =  Timei2; 
load  ran45.mat; 
u4  { : , 5)  =  Timeil; 
y4  ( : , 5)  =  Timei2; 
load  raui46.mat; 
u4(:,6)  =  Timeil; 
y4( : , 6)  =  Timei2; 
load  ran49.mat; 
u4 ( : , 7)  =  Timeil; 
y4{:,7)  =  Timei2; 

%  input  #5 

load  ranSl.mat; 
u5  ( : , 1)  =  Timeil; 
yS  ( : , 1)  =  Timei2; 
load  ran52.mat; 
u5{ : ,2)  =  Timeil; 
y5(:,2)  =  Tiroei2; 
load  r8m53.mat; 
u5(:,3)  =  Timeil; 
y5(:,3)  =  Timei2; 
load  raui54.mat; 
u5 ( : , 4)  =  Timeil; 
y5(:,4)  *  Timei2; 
load  remSS.mat; 
u5(:,5)  =  Timeil; 
y5( : ,5)  =  Timei2; 
load  ran56.mat; 
u5( : , 6)  »  Timeil; 
y5  ( : , 6)  =  Timei2; 
load  ran59.mat; 
u5 ( : , 7)  »  Timeil; 
y5{ : ,7)  «  Timei2; 

%  input  #6 

load  ram6l.mat; 
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u6 ( : , 1)  s  Tlmeil; 
y6  ( : ,  1)  =  Tiinei2; 
load  reua62.n)at; 
u6( : ,2)  =  Timeil; 
y6(:,2)  =  Timei2; 
load  ran63.mat; 
u6(:,3)  =  Timeil; 
y6  ( : , 3 )  =  Timei2 ; 
load  ran64.inat; 
u6 ( : ,4)  =  Timeil; 
y6  ( : , 4)  =  Timei2; 
load  ran65.mat; 
u6 ( : , 5)  =  Timeil; 
y6  ( : , 5)  =  Timei2; 
load  ran66.mat; 
u6 ( : , 6)  =  Timeil; 
y6{:,6)  =  Timei2; 
load  rem69.mat; 
u6 ( : , 7)  =  Timeil; 
y6  ( ; , 7 )  =  Timei2 ; 

%  input  #9 

load  ran91.mat; 
u9  { : , 1)  =  Timeil; 
y9(:,l)  =  Timei2; 
load  ran92.mat; 
u9  ( : , 2 )  =  Timeil ; 
y9  ( : , 2 )  =  Timei2 ; 
load  ran93.mat; 
u9 { : , 3)  =  Timeil; 
y9  ( : , 3 )  =  Timei2 ; 
load  ran94.mat; 
u9 ( : , 4)  =  Timeil; 
y9(:,4)  =  Timei2; 
load  rcui95.mat; 
u9  { : ,5)  =  Timeil; 
y9{:,5)  =  Timei2; 
load  ran96.mat; 
u9(:,6)  =  Timeil; 
y9(:,6)  =  Timei2; 
load  reui99.mat; 
u9(:,7)  =  Timeil; 
y9(:,7)  =  Timei2; 

%  input  #10 

load  ranlOl.mat; 
ulO ( : , 1)  ^  Timeil; 
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ylO(:,l)  *  Timei2; 
load  reml02.inat; 
ul0(:,2)  =  Timeil; 
ylO { : , 2)  =  Timei2 ; 
load  rsuilOS .mat; 
ulO ( : , 3 )  =  Timeil ; 
ylO ( : , 3)  =  Timei2; 
load  ranl04.mat; 
ul0(:,4)  =  Timeil; 
ylO { : , 4)  =  Timei2 ; 
load  reuilOB.mat; 
ul0{:,5)  =  Timeil; 
yl0(:,5)  =  Timei2; 
load  reml06.mat; 
ul0(:,6)  =  Timeil; 
ylO ( : , 6)  =  Timei2 ; 
load  ranl09.mat; 
ul0(:,7)  =  Timeil; 
yl0(:,7)  =  Timei2; 

%  remove  acclercmieter  bias  form  each  measurment 
[r,cl  =  size(yl); 
for  1  =  l:r, 
for  k  =  1:7, 

yl{l,k)  =  yl{l,k)  -  biasCk); 

y2(l,k)  =  y2(l,k)  -  bias{k); 

y3(l,k)  =  y3(l,k)  -  bias(k); 

y4(l,k)  =  y4(l,k)  -  bias{k); 

y5(l,k)  =  y5(l,k)  -  bias(k); 

y6{l,k)  =  y6{l,k)  -  bias(k); 

y9(l,k)  =  y9(l,k)  -  bias(k); 

ylO(l,k)  =  ylO(l,k)  -  bias(k); 
end; 
end; 

%  combine  the  y  and  u  values  to  get  a  mimo  system 
Y  =  yl+y2+y3+y4+y5+y6+y9+yl0; 

U1  =  [ul(:,l)  u2(:,l)  u3(:,l)  u4(:,l)  u5(:,l)  u6(:,l) 
u9(:,l)  U10(:,1)]; 

U2  =  [ul{:,2)  u2(:,2)  u3{:,2)  u4{:,2)  u5(:,2)  u6(:,2) 
u9(:,2)  U10(:,2)]; 

U3  =  [ul(:,3)  u2{:,3)  u3(:,3)  U4(:,3)  u5(:,3)  u6(:,3) 
u9(:,3)  U10{:,3)]; 

U4  =  [ul(:,4)  U2(:,4)  u3(;,4)  u4(:,4)  u5(:,4)  u6(:,4) 
u9(:,4)  U10{:,4)1; 

U5  =  [ul(:,5)  u2(:,5)  u3{:,5)  U4(:,5)  u5(:,5)  u6(:,5) 
u9(:,5)  U10(;,5)]; 

% 


U6  =  [ul(:,6)  U2(:,6)  u3(:,6)  U4(:,6)  u5(:,6)  u6(:,6) 
u9(:,6)  ul0(:,6)]; 

U9  =  [Ul(:,7)  U2(:,7)  u3(:,7)  U4(:,7)  u5{:,7)  u6{:,7) 
u9(:,7)  ul0(:,7)]; 

U  =  U1+U2+U3+U4+U5+U6+U9; 

%  now  begin  system  identification 
p  =  10; 

% [fdl, zmi, shapesl, partial, EMACl]  = 
niimo(Ul,y{ : ,  1)  ,p,  Psan5)le)  ; 

%  [fd2,  zin2,  shapes2,part£a2,EM7lC2]  = 
mimo(U2,  Y( :  ,2)  ,p,  Fsample)  ; 

%  [£d3,  zm3,  shapes3,part£a3,EM7^C3]  = 
mimo{U3,  Y( ; ,  3)  ,p,  Fsample)  ; 

%  [£d4,  zm4,  shapes4,part£a4,EM7^C4]  = 
mimo(U4,  Y( :  ,4)  ,p,  Fsaiiple)  ; 

% [ids, zm5, shapesS, part £a5,EMAC51  = 
mimo(U5, Y{ : , 5) ,p, Fsample) ; 

%  [£d6,  zm6,  shapes6,part£a6,EM;^C€]  = 
mimo{U6,  Y{ : ,  6)  ,p,  Fsanple) ; 

%  [£d9,  zm9,  shapes9,part£a9,  EM71C9]  = 
mimo(U9,  Y( : ,  7)  ,p,  Fsanple)  ; 

[£d,  zm, poles,  shapes, partfac]  =  mdLmo(U,Y,p,  Fsanple)  ; 
cut  =  0.001; 

%  compute  natural  frequency 

%  evaluate  model 

wdm=2*pi*£d; 

for  1  =  1 ; length ( fd) , 

wnrnd)  =wdm(l)  /sqrt  (l-zm(l)  *zm{l) ) ; 
end; 

%  Create  state  space  model  A  matrix 
k  =  1; 

for  1  =  1 : length (wnm) , 

z  (k)  =exp(  ( -zm(l)  *wnm(l)  +j*wdm{l) )  /Fsanple)  ; 
z (k+l) =conj (z (k) ) ; 

B(k, : ) =partfac{l,  ; ) ; 

B(k+1, : ) =conj {B{k,  : ) )  ; 
cm( : ,  k)  =shapes  ( : ,  1)  ; 
cm( :  ,k+l)  =conj  (cm( ;  ,k) )  ; 
k  =  k  +  2; 
end; 

am  =  diag(z)  ; 

[rl,cl]  =  size(am); 

[atenp,bm]  =  c2d(am,  B,  l/Fsanple) ; 
dm  =  zeros (7, 8) ; 

%  transform  to  real  matrices 


97 


T1  =  [1  l;j  -j]; 

T  =  Tl; 

for  k  =  4:2:rl, 

[r,c]  =  size(T); 

T  =  [T  zeros (r, 2) ; zeros (2, c)  Tl] ; 
end; 

Tinv  =  inv(T) ; 

A  =  T*am*Tinv; 

B  =  T*bm; 

C  =  cm*Tinv; 

D  =  dm; 

save  tmabc.mat  A  B  C  D 

save  tmresul . mat  wnm  fd  zm  shapes  partfac 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  tmdl 

% 

%  This  file  creates  the  mcxlel  C 
%  matrix  for  the  PACOSS.  The  A 

%  euid  B  matrices  are  confuted  from 

%  either  the  ERA  or  BLS  methods. 

%  This  script  file  uses  PKM  to 

%  approximate  the  C  matrix  one  row 

%  at  a  time.  It  only  uses  the  data  from 

%  the  colocated  actuator/accelerometer 

%  pairs. 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  input  the  PACOSS  data  needed  to  generate 

%  the  C  matrix 

%  Input  #1 

load  trcinll.mat; 

ul(:,l)  =  Timei4 (45:512)  ; 

yl(:,l)  =  Timeil(45 : 512)  ; 

%  input  #2 

load  tran22.mat; 

U2(:,2)  =  Timei4 (45 ; 512)  ; 
y2(:,2)  =  Timeil(45:512)  ; 

%  input  #3 

load  tran33.mat; 
u3(:,3)  =r  Timei4  (45:512)  ; 
y3(:,3)  =  Timeil (45 : 512)  ; 

%  input  #4 

load  tran44 . mat ; 
u4(:,4)  =  Timei4(45:512) ; 
y4(:,4)  =  Timeil (45 : 512) ; 

%  input  #5 

load  tran55.mat 
u5(:,5)  =  Timei4(45:512) ; 
y5(:,5)  =  Timeil (45 : 512)  ; 

%  input  #6 

load  trams 6. mat; 

U6(:,6)  =  Timei4 (45:512)  ; 
y6(:,6)  =  Timeil (45 : 512)  ; 

%  input  #9 

load  tran99.mat; 


99 


u9(:,7)  =  Tiinei4(45:512)  ; 
y9{:,7)  =  Timeil(45:512)  ; 

%  input  #10 

load  tremlOlO .mat; 

U10(:,8)  =  Timei4{45:512) ; 
yl0{:,8)  =  Timeil(45:512) ; 

%  subtract  out  the  accelerometer  bias 

%  the  mean  Is  used  because  the  actual  system 

%  has  no  net  motion,  so  the  mean  is  zero 

yl(:,l)  =  yl(:,l)  -  mean(yl  ( : ,  1) )  ; 

y2(:,2)  =  y2{:,2)  -  mean(y2  ( :  ,2) )  ; 

y3(:,3)  =  y3(:,3)  -  mean(y3  ( : ,  3) )  ; 

y4(:,4)  =  y4{;,4)  -  meeUl  (y4  ( : ,  4) )  ; 

y5(:,5)  =  y5{:,5)  -  mean(y5  ( : ,  5) )  ; 

y6(:,6)  =  y6(:,6)  -  mean(y6  ( : ,  6) )  ; 

y9(:,7)  =  y9(:,7)  -  mecui(y9  ( : ,  7) )  ; 

yl0{:,8)  =  yl0(:,8)  -  meeui(yl0  ( : ,  8) )  ; 

%  use  PEM  to  estimate  C  matrix 
%  load  the  A,  B,  euid  C  matrices 
load  tabc.mat 
%load  tmabc.mat 
T  =  1/25.6; 
frowa,cola]  =  size (A); 

[rowb,colb]  =  size(B); 

[rowc,colc]  =  size(C); 

ai  =  A; 

bil  =  B(:,l); 

bi2  =  B{: ,2) ; 

bi3  =  B(:,3); 

bi4  =  B{: ,4) ; 

bis  =  B(; ,5) ; 

bi6  =  B(:,6) ; 

bi9  =  B{:  ,7)  ; 

bilO  =  B(: ,8) ; 

ci  =  nan*ones (1, rowa) ; 

di  =  zeros ( 1, 1); 

ki  =  zeros (rowa, 1) ; 

xOi  =  zeros (rowa, 1) ; 

msil  =  modstznac(ai,bil,ci,di,ki,xOi) ; 

msi2  =  modstruc(ai,bi2, ci,di,ki,xOi) ; 

msi3  =  modstruc (ai,bi3, ci,di,ki,xOi)  ; 

msi4  =  modstruc(ai,bi4, ci,di,ki,xOi)  ; 

msi5  =  modstruc (ai, bis, ci,di,ki,xOi) ; 

msi6  =  modstruc(ai,bi6, ci,di,ki,xOi) ; 

msi9  =  modstruc(ai,bi9, ci,di,ki,xOi) ; 
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msllO  «  inod8truc(al,bllO,cl,dl,ki,xOl)  ; 

%  initial  guess  for  the  C  matrix 

paxval  s  C ( 1 , : ) ; 

parva2  =  C(2,:); 

parva3  =  C{3,:); 

parva4  =  C{4,:); 

parvaB  *  C ( 5 , : ) ; 

parvae  *  C ( 6 , : ) ; 

parva9  =  C ( 7 , : ) ; 

parvalO  =  C ( 8 , : ) ; 

lambdi  »  [] ; 

thil  =  ms2th(insil, ’d*  ,parvai,  lambdi, T)  ; 
thi2  =  m82th(m8i2, 'd' ,parva2, lambdi, T) ; 
thi3  =  ms2th(msi3, 'd' ,parva3,lainbdi,T)  ; 
thi4  =  ms2th(msi4, 'd' ,parva4, lambdi, T) ; 
this  =  ms2th(m8i5, 'd* , parvaS , lambdi , T) ; 
thi6  =  ms2th(msi€, 'd' ,parva6, lambdi, T) ; 
this  =  ms2th(m8i9, 'd' ,parva9, lambdi, T) ; 
thilO  =  ms2th(msil0, 'd' , parvalO, laxnbdi,T) ; 
index  =  [1 :  length  (peunral)  ] ; 

%  perform  system  ID 

thl  =  pem(tyl(:,l)  ul (:, 1) ], thil, index, -l, le- 10, -l, -l,T) 
th2  =  pem([y2(:,2)  u2 ( : , 2) ] , thi2, index, -1, le-10, -1, -1,T) 
th3  =  pem([y3(:,3)  u3 3) ], thi3, index, -1, le-10, -1, -1, T) 
th4  =  pem(ty4(:,4)  u4( : ,4) ] ,thi4, index, -1, le-10, -l, -l,T) 
ths  =  pem(ty5(:,5)  u5 (:, 5) ], this, index, -1, le-10, -1, -1,T) 
the  =  pem(ty6(;,6)  u6 (:, 6) ], thie, index, -l, le-10, -l, -l,T) 
th9  =  pem(ty9(:,7)  u9 (:, 7) ], this, index, -1, le-10, -1, -1,T) 
thlO  =  pem( [ylO ( : , 8)  ulO {:, 8) 1 , thilO, index, -1, le-10, -1, - 
1,T)  ; 

%  convert  theta  to  SS  format 
[am,bnil, cml,dm, kml,x011  «  th288(thl); 

[am, bin2, cm2,dm, kml,x0l]  >  th288(th2); 

[am,bm3, cm3,dm,kml,x01]  «  th288(th3); 

[am» bm4,  cm4,dm,kinl,x01]  »  th288(th4); 

[am,bm5,cm5,dm, kml,x0l]  »  th288(th5); 

[am, bme, cm6,dm, kinl,x0l]  «  th288(th6); 

[am,bni9, cm9,dm, kml,x01]  *  th288(th9); 
[am,binl0,cml0,dm,kinl,x01]  >  th2ss(thl0)  ; 

%  ccnibine  c  matrices  to  get  MIMO  C  matrix 

am  «  A; 

bm  «  B; 

cm  «  [cml;cm2;cm3;cm4;cm5;cm€;cxn9;cml0]  ; 


dm  =  D; 

save  tmodl.mat  am  bm  cm  dm 
%save  tmnodl.mat  am  bm  cm  dm 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  tfcoipi 

% 

%  This  file  conpares  the  TF  of  the  truth  system 
%  and  the  model  system.  Magnitues  plots  are  in 
%  db,  and  phase  plots  are  in  degrees.  The  TF's 
%  are  from  l  to  10  Hz 

% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% 

%  load  the  PACOSS  TF  data  and 
%  generate  TF  plots  from  PACOSS  data 
FreqV  =  0:0.025:10; 

%  input  #1 

load  tfvl 

TFll  =  20*logl0(abs(vll(41:40l) ) ) ; 

TF12  =  20*logl0(abs(vl2(41:40l))) ; 

TF13  =  20*logl0 {abs(vl3 (41:401) )) ; 

TF14  =  20*logl0 (abs (Vl4 (41:401) )) ; 

TF15  =  20*logl0 (abs(vl5 (41:401) )) ; 

TP16  =  20*logl0(abs(vl6(41;40l) ) ) ; 

TF19  =  20*logl0 (abs (vl9 (41:401) )) ; 
phasell  = 

(180/pi) *unwrap(atan2 (imag(vll (41:401) ) ,real (vll (41:401) ) ) ) 
phasel2  = 

(180/pi)  *unwrap(atcui2  (imag(vl2  (41:401) ) ,  real  (vl2  (41:401) )) ) 
phasel3  = 

(180/pi)  *unwrap(ateui2  (imag(vl3  (41:401) )  ,real  (vl3  (41:401) ) ) ) 
phasei4  = 

(180/pi) *unwrap(atan2 (imag(vl4 (41:401) ) ,real (vl4 (41:401) ) ) ) 
phaselS  = 

(l8C/pi) *unwrap(atan2 (imag(vl5(4l:40l) ) ,real (vl5(4l:40l) ) ) ) 
phasei6  = 

(180/pi)  *unwrap(atcui2  (imag(vl6  (41:401) )  ,real(vl6(4l:40l) ) ) ) 
phasei9  = 

(180/pi) *unwrap(atan2 (imag(vl9 (41:401) ) ,real (vl9 (41:401) ) ) ) 

%  input  #2 

load  tfv2 

TP21  =  20*logl0 (abs(v21(41:401) ) ) ; 

TF22  =  20*logl0  (cJas  (v22  (41:401) ))  ; 

TP23  =  20*logl0(abs(v23(41:401) ) ) ; 

TF24  *  20*logl0 (abs (v24 (41:401) )) ; 
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TP25  =  20*logl0 (abs(v25 (41:401) )) ; 

TP26  =  20*logl0(abs(v26(41:40l) ) ) ; 

TP29  =  20*logl0 (abs{v29 (41:401) )) ; 
phase2l  s: 

(180/pi)  *unwrap(atem2  (linag(v21  (41:401) )  ,real  (v2l  (41:401) )))  ; 
phase22  s 

(180/pi)  *unwrap(atan2  (iniag(v22  (41:401) ) ,  real  (v22  (41:401) ) ) ) ; 
phase23  = 

(180/pi) *unwrap(atan2 (imag(v23 (41:401) ) , real (v23 (41:401) ) ) ) 
phase24  = 

(180/pi)  *unwrap(atcui2  (iinag(v24  (41:401) )  ,  real  (v24  (41:401) ) ) )  ; 
pha8e25  = 

(180/pi)  *unwrap(atoUx2  (iinag(v25  (41:401) ) ,  real  (v25  (41:401) ) ) )  ; 
phase26  = 

(180/pi)  *tinwrap(atan2  (iinag(v26 (41:401) )  ,real  (v26  (41:401) ) ) )  ; 
phase29  = 

(180/pi)  *unwrap(ateu:i2  (imag(v29  (41:401) )  ,real  (v29  (41:401) ) ) )  ; 

%  input  #3 

load  tfv3 

TP31  =  20*logl0(abs(v31(41:401) ) ) ; 

TP32  =  20*logl0(abs(v32(41:401) ) ) ; 

TP33  =  20*logl0 (abs(v33 (41:401) )) ; 

TP34  =  20*logl0(abs(v34(41:401) ) ) ; 

TP35  =  20*logl0 (abs(v35 (41:401) )) ; 

TP36  =  20*logl0(ab8(v36(41:401) ) ) ; 

TP39  =  20*logl0(ab8(v39(41:401) ) ) ; 
pha8e31  - 

(180/pi)  *unwrap(atan2  (iinag(v3l(41:401) )  ,real  (v3l(4l:401) ) ) ) ; 
pha8e32  = 

(180/pi)  *unwrap(ataui2  (iinag(v32  (41:401) )  ,real  (v32  (41:401) ) ) )  ; 
pha8e33  - 

(180/pi)  *unwrap(ataui2 (invag(v33 (41:401) )  ,real (v33(41:401) ) ) )  ; 
pha8e34  = 

(180/pi)  *unwrap(atan2  (iinag(v34  (41:401) )  ,  real  (v34  (41:401) ) ) )  ; 
pha8e35  = 

(180/pi)  *unwrap(atcui2  (iinag(v35  (41:401) ) ,  real  (v35  (41:401) ) ) )  ; 
pha8e36  = 

(180/pi)  *unwrap(atan2  (iniag(v36  (41:401) )  .real  (v36  (41:401) ) ) )  ; 
pha8e39  s 

(180/pi)  *unwrap(atan2  (iinag(v39  (41:401) ) ,  real  (v39  (41:401) ) ) )  ; 
%  input  #4 

load  tfv4 

TP41  =  20*logl0 (ab8 (v41 (41:401) )) ; 

TP42  »  20*logl0 (ab8(v42 (41:401))); 
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TF43  =  20*logl0(abs(v43(41:40l))); 

TF44  »  20*logl0 (ab8(v44 (41:401) )) ; 

TF45  =  20*logl0 (abs(v45 (41:401) >) ; 

TF46  =  20*logl0(ab8(v46(41:401) ) ) ; 

TF49  =  20*logl0(ab8(v49(41:401) ) ) ; 
pha8e41  = 

(180/pi)  *unwrap(ataui2  (iinag(v4l(4l:40l) )  ,real  (v4l  (41:401) ) ) )  ; 
pha8e42  - 

(180/pi)  ♦unwrap  ( atan2  (iinag(v42  (41:401) ) ,  real  (v42  (41:401) ) ) )  ; 
pha8e43  = 

(l80/pi)  *unwrap(ataui2  (ixnag(v43  (41:401) )  ,real(v43(41:401) ) ) )  ; 
pha8e44  = 

(180/pi)  ♦iin%in:ap(atan2  (imag(v44  (41:401) ) ,  real  (v44  (41:401) ) ) )  ; 
phase45  » 

(180/pi)  *unwrap(atan2  (iinag(v45  (41:401) ) ,  real  (v45  (41:401) ) ) )  ; 
pha8e46  s 

(180/pi)  ♦unwrap  (atan2  (iinag(v46  (41:401) ) ,  real  (v46  (41:401) ) ) )  ; 
pha8e49  « 

(180/pi)  ♦uniirrap(atan2  (iniag(v49  (41:401) ) ,  real  (v49  (41:401) ) ) )  ; 
%  input  #5 

load  tfv5 

TF51  =  20^1ogl0 (ab8(v51 (41:401) )) ; 

TFS2  =  20^1ogl0(ab8(v52(41:401) ) ) ; 

TF53  »  20^1ogl0(ab8(v53(41:401) ) ) ; 

TF54  =  20^1ogl0(ab8(v54(41:401) ) ) ; 

TF55  =  20^1ogl0(ab8(v55(41:401)) ) ; 

TF56  =  20^1ogl0(ab8(v56(41:401) ) ) ; 

TF59  =  20^1ogl0 (ab8(v59 (41:401) )) ; 
phaeeSl  = 

(180/pi)  ♦unwrap (atan2  (iinag(v5l(4l:40l) )  ,real(v5l(4l:401) ) ) )  ; 
pha8e52  = 

(180/pi)  ♦unwrap ( ataui2  (iniag(v52  (41:401) )  ,  real  (v52  (41 : 401) ) ) )  ; 
pha8e53  = 

(180/pl)  ♦unwrap  ( atan2  (iinag(v53  (41:401) ) ,  real  (v53  (41:401) ) ) )  ; 
pha8e54  - 

(180/pi)  ♦unwrap  (atcm2  (iniag(v54  (41:401) ) ,  real  (v54  (41:401) ) ) )  ; 
phaaeSS  > 

(180/pi)  ♦unwrap (ataui2  (iinag(v55  (41:401) )  ,real(v55(41:401) ) ) )  ; 
pha8e56  « 

(180/pi)  ♦\inwrap(atan2  (iinag(v56  (41:401) ) ,  real  (v56  (4l :  401) ) ) )  ; 
pha8e59  « 

(180/pi)  ♦unwrap ( at am2  (iinag(v59  (41:401) )  ,real  (v59  (41:401) ) ) )  ; 
%  input  #6 

load  tfv6 

TF61  -  20^1ogl0 (ab8(v61 (41:401))); 
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TP62  =  20*logl0{abs(v62{41;401))); 

TP63  =  20*logl0(ab8(v63(41:401) ) ) ; 

TP64  =  20*logl0 (abs (v64 (41:401) )) ; 

TP65  *  20*logl0(abs(v65(41:401) ) ) ; 

TP66  =  20*logl0(ab8(v66(41:40l) ) ) ; 

TP69  =  20*logl0 (ab8(v69 (41:401) )) ; 
pha8e6l  - 

(180/pl)  *iinwrap(ataui2  (linag(v6l  (41:401) )  ,real  (v6l  (41:401) ) ) ) ; 
pha8ee2  = 

(180/pi)  *unwrap(atoUi2  (linag(v62  (41:401) )  ,real  (v62(41:401) ) ) ) ; 
pha8e63  » 

(180/pi)  *unwrap ( atan2  (iinag(v63  (41:401) ) ,  real  (v63  (41:401) ) ) ) ; 
pha8e€4  s 

(180/pi)  *unwrap(atan2  (iinag(v64  (41:401) ) ,  real  (v64  (41:401) ) ) )  ; 
phaaeCS  » 

(180/pi)  *iin«rrap(atan2  (imag(v65  (41:401) ) ,  real  (v65  (41:401) ) ) ) ; 
pha8e€6  « 

(180/pl)  *xinwrap(atan2  (iinag(v66  (41:401) )  ,  real  (v66  (41:401) ) ) )  ; 
pha8e69  « 

(180/pi)  *uinwrap(atam2  (iinag(v69  (41:401) ) ,  real  (v69  (41 : 401) ) ) )  ; 

%  input  #9 

load  tfv9 

TP91  =  20*logl0(ab8(v91(41:401))) ; 

TF92  =  20*logl0(ab8(v92(41:401) ) ) ; 

TP93  =  20*logl0(ab8(v93(41:401) ) ) ; 

TP94  =  20*logl0 (ab8(v94 (41:401) )) ; 

TP95  =  20*logl0(ab8(v95(41:401) ) ) ; 

TP96  =  20*logl0 (ab8(v96 (41:401) )) ; 

TP99  =  20*logl0 (ab8(v99 (41:401) )) ; 
pha8e9l  « 

(180/pi)  *un«nrap(atan2  (ln)ag(v91  (41:401) ) ,  real  (v91  (41:401) ) ) )  ; 
pha8e92  « 

(180/pi)  *\inwrap ( atan2  (iinag(v92  (41:401) )  ,real  (v92  (41:401) ) ) )  ; 
pha8e93  « 

(180/pi)  *iinwrap(atan2  (iinag(v93  (41:401) )  ,real  (v93  (41:401) ) ) )  ; 
pha8e94  « 

(180/pi)  *un«rrap(atan2  (iinag(v94  (41:401) )  ,real  (v94  (41:401) ) ) )  ; 
phase95  « 

(180/pl)  *unwrap(atan2  (linag(v95  (41:401) ) ,  real  (v95  (41:401) ) ) )  ; 
pha8e96  « 

(180/pi)  *unwrap(atan2  (iinag(v96  (41:401) )  ,real(v96(41:401) ) ) )  ; 

pha8e99  « 

(180/pi)  *unwrap(ataua2  (iinag(v99  (41: 401) )  ,real(v99(4l:40l) ) ) )  ; 
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%  input  #10 

load  tfvio 

TPlOl  =  20*logl0 (ab8{vl01 (41:401) )) ; 

TP102  *  20*logl0 (ab8(vl02 (41:401) )) ; 

TP103  =  20*logl0(ab8(vl03(41:401) ) ) ; 

TP104  =  20*logi0 (ab8(vl04 (41:401) )) ; 

TP105  =  20*logl0(ab8(vl05(41:401) ) ) ; 

TP106  =  20*logl0(ab8(vl06(41:401) ) ) ; 

TP109  =  20*logl0 (ab8 (vl09 (41:401) )) ; 

TPlOlO  =  20*logl0 (ab8(vlOlO (41:401) )) ; 
pha8el01s  (180/pi)  *uniin:ap(atan2  (iinag(vlOl  (41:401) ) ,  real  (vlOl  ( 
41:401) ))); 

pha8el02s  (180/pl)  *unwrap(atan2  (linag(vl02  (41:401) ) ,  real  (vl02  ( 
41:401) ) ) )  ; 

phaeelOls  (180/pl)  ♦un%frap(atan2  (iinag(vl03  (41:401) ) ,  real  (vi03  ( 
41:401) ) ) )  ; 

pha8el04s (180/pi)  *unwrap(atan2 (iinag(vl04 (41:401) )  ,real(vl04( 
41:401) ) ) )  ; 

phaeelOSs  (180/pi)  ♦unwrap  ( ataui2  (iinag(vl05  (41:401) )  ,real  (vl05  ( 
41:401) ) ) )  ; 

pha8el06s  (180/pi)  *un%n:ap(atan2  (iinag(vl06  (41:401) ) ,  real  (vi06  ( 
41:401) ) ) )  ; 

pha8el09=  (180/pi)  ♦unwrap (atan2  (iioag(vl09  (41:401) )  ,real  (vl09  ( 
41:401) ) ) )  ; 

pha8el010>:  (180/pl)  ♦\inwrap(atan2  (iinag(vl010  (41:401) ) ,  real  (vlO 
10(41:401) ) ) ) ; 

%  load  model  A,  B,  C,  emd  D  roatricee 
load  tmodl 
%load  tmmodl 

%  generate  bode  plot  of  eyatem  model 
wrad  =  2^pi^PregV; 

%  reduce  gain  of  eyatem 
cm  *  O.Ol^cm; 

[magml,  phasl]  «dbode  (am, bm,  cm, dm,  1/51 .2,1,  wrad)  ; 

[magm2,pha82]  sdbode(am,bin,cm,dm,  l/51.2,2,wrad)  ; 

[magm3 ,  pha83 ]  -dbode  ( am,  bm,  cm,  dm,  l/5l .  2 , 3 ,  wrad)  ; 

[magm4,pha84l «dbode(am,bm, cm, dm,  1/51.2, 4, wrad)  ; 

[magm5,pha85]  ^dbode ( am, bm, cm, dm,  1/51.2, 5, wrad) ; 

[magm6 , phaaG]  «dbode (am, bm,  cm, dm,  l/5l .2,6,  wrad) ; 

[magm9,pha89]  «dbode(am,bm,  cm, dm,  1/51.2, 7, wrad) ; 
[ffiaginl0,pha8l0]  *dbode(am,bm,cm,dm,  1/51.2, 8,  wrad) ; 

%  compare  truth  suid  model 

plot  (Pre<JV(41: 401)  ,  [TPll,20^1ogl0(magmi(4l:401, 1) )  1 ) 
aave  t final. mat  aun  bm  cm  dm 
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